## Description

1 Revisiting Image Convolution & Deconvolution

Let’s revisit and generalize the image convolution and deconvolution problems we saw in Assignment 1. This

time, instead of only consider a discrete 1D horizontal blur kernel, we will consider more general 2D blurs

(i.e., horizontal and vertical).

To illustrate the forward blurring process we extend our previous diagrams, however the “sweep and stamp”

process still applies. Our diagram below illustrates a few instances of applying a blur kernel (i.e., with 9

elements, which we’ve numbered for clarity: ❶, ❷, , ❾) to a input image, and generating a

output image, e.g.,

⬚ ⬚ ⬚ ⬚ ⬚

❶❷❸ ⬚ ⬚ ⬚ ⬚ ⬚

❹❺❻ ⬚ ⬚ ⬚ ⬚

❼❽❾ ⬚ ⬚ ⬚ ⬚ ⬚

⬚ ⬚ ⬚ ⬚ ⬚

3 × 3

… 5 × 5 5 × 5

You can assume that we will only test square kernels, input/output image resolutions greater than

or equal to the kernel size, and only square images. That being said, writing your code in a manner

that doesn’t assume square kernels can help with debugging/testing.

ꋏ

The kernel is centered — this time horizontally and vertically — on an input pixel at the location of the

corresponding output pixel (illustrated as a circle with black contour and grey fill), and if portions of the

kernel happen to fall outside the bounds of the input image, they wrap around (vertically and horizontally,

now).

By “sliding” the kernel along each pixel of the uncorrupted input image, weighting the underlying

source/input pixels by the overlapping kernel values and summing across these weighted input pixels, we

construct the final blurred image one pixel at a time:

❺❻ ❹ ⬚ ⬚ ⬚ ⬚ ❹❺❻ ⬚ ⬚ ⬚ ⬚ ❹❺❻ ⬚ ⬚ ⬚ ⬚

❽❾ ❼ ⬚ ⬚ ⬚ ⬚ ⬚ ❼❽❾ ⬚ ⬚ ⬚ ⬚ ⬚ ❼❽❾ ⬚ ⬚ ⬚ ⬚ ⬚

⬚ ⬚ ⬚ ⬚ ⬚ , ⬚ ⬚ ⬚ ⬚ ⬚ , ⬚ ⬚ ⬚ ⬚ ⬚ , ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚

❷❸ ❶ ⬚ ⬚ ⬚ ⬚ ⬚ ❶❷❸ ⬚ ⬚ ⬚ ⬚ ⬚ ❶❷❸ ⬚ ⬚ ⬚ ⬚ ⬚

❻ ❹❺ ⬚ ⬚ ⬚ ⬚ ❾ ❼❽ ⬚ ⬚ ⬚ ⬚ ⬚

❾ ❼❽ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ . . . , ⬚ ⬚ ⬚ ⬚ ⬚ , . . . , ⬚ ⬚ ⬚ ⬚ ⬚

⬚ ⬚ ⬚ ⬚ ⬚ ❸ ❶❷ ⬚ ⬚ ⬚ ⬚ ⬚

❸ ❶❷ ⬚ ⬚ ⬚ ⬚ ⬚ ❻ ❹❺ ⬚ ⬚ ⬚ ⬚

Recall, as in Assignment 1, while the diagrams above illustrate the blurring process in 2D and for 2D images

and kernels, you will express this linear map as a matrix — and thus — must treat the input (and output)

images as “flattened” 1D vectors. Treating this change of indexing from 2D pixel coordinates to flattened 1D

“image” vector elements is one part of the learning objective for this task.

1.1 Isotropic Gaussian Blur Kernel

One common blur kernel is the discrete isotropic Gaussian kernel. For one-dimensional blurs, the continuous

kernel is defined as

where is the standard deviation of the Gaussian, which defines its fall-off rate, and the continuous kernel is

defined over the entire real line, i.e., for . In practice, one way to take this continuous kernel

and discretize it to have length is to evaluate the kernel at integer lattice values .

For two-dimensional blur kernels, the continuous isotropic Gaussian is

Deliverable 1 [25%]

Construct the blur matrix: Complete the implementation of CreateBlurMatrix. It takes as input

a numpy.array with the 2D blur kernel elements, as well as two positive integers denoting

the width and height of the image the matrix will be applied to — you can assume that the

corrupted and uncorrupted images have the same dimensions, and that all images are square (i.e.,

width = height). The routine returns the blur matrix as a numpy.array with shape (width

height, width height).

☑

(푘, 푘)

×

×

퐺1D (푥; 휎) = 1

√‾2‾휋‾휎

푒

− 푥2

2휎 2 (1)

휎

−∞ ≤ 푥 ≤ ∞

푘 푥 ∈ {−⌊ ⌋, … , ⌊ ⌋} 푘

2

푘

2

and its discretized version is obtained by evaluating on a 2D integer lattice with

and . Note that can be decomposed as the

product of two 1D Gaussians, a fact you can optionally leverage to simplify your implementation.

1.2 Conditioning & Regularization

As in Assignment 1, we provide you with a benchmark image set. Here, in addition to providing the clean and

blurred images, we also provide a version of the blurred image polluted with subtle noise. We will use this

noisy blurred image to explore the conditioning of your blur operator.

In the current scenario, with a fully-constrained blur matrix and blurred image solving for the leastsquares unpolluted image as

amounts to solving the original linear system as . For many blurring kernels, the blur matrix can

become ill-conditioned, meaning the solution can be catastrophically sensitive to minor deviations in the

input .

We can employ regularization strategies to try to mitigate the poor conditioning of the system. One common

strategy is Tikhonov regularization, which aims instead to solve the modified least-squares problem in

where is the Tikhonov regularization parameter that allows you to control the degree of regularization. This

regularized least-squares problem can be expressed as the solution to the following augmented system of

linear equations:

퐺2D (푥, 푦; 휎) = 1

2휋휎2

푒

− 푥2+푦2

2휎 2 (2)

푘 × 푘 (푥, 푦)

푥 ∈ {−⌊푘/2⌋, … , ⌊푘/2⌋} 푦 ∈ {−⌊푘/2⌋, … , ⌊푘/2⌋} 퐺2D

Deliverable 2 [10%]

2D Gaussian Blur Kernel: Complete the implementation of Gaussian2D to generate the 2D

discretized Gaussian blur kernel. It takes as input the (odd) kernel extent (i.e., for a square 2D

kernel of size ) and the standard deviation for the isotropic Gaussian sigma. It outputs a

numpy.array of the 2D Gaussian evaluated at integer lattice locations.

☑

푘 푘

(푘, 푘) 휎

(푘, 푘)

Consider different ways that you can test your Gaussian2D implementation, as well as how it can

be used to test your implementation of CreateBlurMatrix.

ꋏ

퐀 퐛

퐱

min ‖퐛 − 퐀퐱 퐱 ‖2

2 (3)

퐱 = 퐀 퐛 −1

퐛

min (‖퐛 − 퐀퐱 + ‖퐱 ) = , 퐱 ‖2

2 휆2 ‖2

2 min

퐱 [ ] − [ ] 퐱

‖

‖

‖ 퐛

0

퐀

휆퐈

‖

‖

‖

2

2

(4)

휆

(퐀 퐀 + 퐈) 퐱 = 퐛 . T 휆2퐈

T 퐀T (5)

Tikhonov regularization is arguably the most basic regularization scheme but, even then, a judicious choice

of its parameter can yield much better results than a more naïve setting. One generalization of the

Tikhonov formulation,

admits an additional operator — a so-called regularization matrix — that can be designed to specialize the

regularization based on the needs of a specific application. In the case of image deblurring, some choices for

that can outperform vanilla Tikhonov regularization include local derivative and/or Laplacian operators.

2 Solving Linear Systems with Singular Value Decomposition

Let’s explore how singular value decompositions (SVDs) can be used to solve (implicitly regularized) linear

systems. As a warm up, we will implement a simple singular vector solver, before moving on to the

deconvolution problem using numpy’s built-in SVD algorithm. Recall that the SVD of is and has

non-zero singular values.

Deliverable 3 [10%]

Regularized deblurring: Complete the implementation of DeblurImage. It takes as input the

unregularized blur matrix, the (potentially noisy) blurred image, and the Tikhonov regularization

constant lambda_. The function should construct the augmented/regularized linear system of

equations and solve them (using whichever numpy routine you like) to solve for and output the

(unflattened, 2D) deblurred image.

☑

휆

• You may find it useful to test your routines on smaller images.

• You may want to complete the rest of the assignment before returning to complete the next

deliverable.

ꋏ

휆

min (‖퐛 − 퐀퐱 + ‖퐃퐱 ) = , 퐱 ‖2

2 휆2 ‖2

2 min

퐱 [ ] − [ ] 퐱

‖

‖

‖ 퐛

0

퐀

휆퐃

‖

‖

‖

2

2

(6)

퐃

퐃

Deliverable 4 [15%]

Deblurring competition: DeblurImage implements the basic Tikhonov regularization scheme,

however other regularization methods for deblurring exist. Complete the implementation of

DeblurImageCompetition to implement any set of regularization approaches with the goal of

generating your best possible deblurred image result — here, you have full latitude in how you go

about regularizing, augmenting, etc. the linear system in order to improve the deblurring result.

The function signature is the same as DeblurImage from A1 (i.e., the same as DeblurImage

above, but without lambda_; you should hardcode any regularization constants you use within

your DeblurImageCompetition function body). Your grade in this deliverable will be based on

how well your deblurring works compared to your colleagues’, and how well you document your

solution. Deblurring error will be computed as the difference between your deblurred image

and the ground truth deblurred (i.e., source) image.

☑

퐿2

퐀 퐀 = 퐔횺퐕T

푘 ≤ min(푚, 푛)

2.1 Power Iteration

We discussed how the Power Iteration Method could be used to compute the eigenvectors associated to the

largest and smallest eigenvalues in a non-singular matrix . We can extend this approach to compute the left

and right singular vectors associated to the largest and smallest singular values of , by applying the method

to two separate eigenvector settings:

by applying the power iteration method to you can compute the left singular vectors in

associated to the largest ( for ) or smallest ( for ) singular values, and

by applying the power iteration method to you can compute the left singular vectors in

associated to the largest ( for ) or smallest ( for ) singular values.

2.2 Rank-reduced Deconvolution

We can use the SVD to solve linear systems. Consider taking the SVD of your blur matrix as and

then solving for in the unregularized deblurring problem. Clearly if we use the full-rank

reconstruction of as we will run into conditioning-related instabilities when, e.g., the blurred image

is polluted by noise.

Instead of using the full rank SVD reconstruction of , you can implement a rank-reduced reconstruction

that retains only a subset (i.e., a percentage up to 100%) of the system’s “energy”. To do so whilst retaining

of the system’s energy, you can form an effective rank-reduced approximation of that only

retains those singular values above the -percentage threshold — and zeros out the remaining singular

values.

Specifically, your rank-reduced approximation of will include only those singular values that

satisfy this inequality:

Your rank-reduced has and can be used to solve

Here, we assume that (i.e., ), otherwise the degree to which is zero-padded would differ.

Keep in mind that, while implementing Equation as-is will not lead to incorrect results, there are several

avenues for optimizing its performance without sacrificing any accuracy — in fact, some optimizations may

improve the numerics of the solution.

퐀

퐀

퐀퐀T 퐔

퐮1 휎1 퐮푘 휎푘

퐀 퐀T 퐕

퐯1 휎1 퐯푘 휎푘

Deliverable 5 [15%]

Power Iteration: Complete the implementation of PowerMiniSVD to return the left and right

singular vectors associated to the largest singular value, and . The function takes an input

matrix and number of power iterations as input.

☑

퐮1 퐯1

퐀 = 퐔횺퐕T

퐱 = 퐀 퐛 −1

퐀 퐔횺퐕T

퐛

퐀

0 ≤ 푅 ≤ 1 퐀

푅

퐀 {휎1, … , 휎푟}

(∑ / ) ≤ 푅 . 푟

푖=1 휎2

푖 ∑푘

푗=1 휎2

푗 (7)

퐀̂= 퐔횺̂

퐕T 횺 = diag( , … , , 0, … , 0) ̂ 휎1 휎푟

퐱 = 퐀 퐛 = 퐕 퐛 . ̂

−1

횺̂

−1

퐔T (8)

푟 < 푘 푅 < 1 횺̂

8

3 You’re Done!

Congratulations, you’ve completed the 3rd assignment. Review the submission procedures and guidelines at

the start of the Assignment 0 handout before submitting your Python script file assignment solution.

formatted by Markdeep 1.13

Deliverable 6 [25%]

SVD Deblurring: Complete the implementation of DeblurImageSVD which takes as input the SVD

factors of a blur matrix (computed outside the function), the blurred (and potentially noisy) input

image, the amount of energy to retain in the rank-reduced reconstruction of the blur matrix, and

returns the deblurred image as in Equation .

☑

푅

8

Some notes:

• Computing the SVD is expensive. You may want to implement a caching scheme in your test

code using numpy.save and numpy.load

• Do your rank-reduced image deblurring results require additional explicit regularization? You

don’t have to actually answer this in your solution, but it’s worthwhile reflecting on.

ꋏ

✒