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.
ꋏ
✒
