CMDA 3606, Spring 2021 Homework 5 Dr. Chung Due Thursday March 18, 2021 at 2:00pm
The goal of this project is to derive the discrete system of linear equations that represents discrete 1D convolution1. Recall that in the continuous formulation, convolution has the form
1 0
where we have changed variable names for notational convenience. In the discrete formula- tion, each value of the function b(s) is essentially a weighted average of the values of x(t), where the weights are given by the kernel function or point spread function (PSF), p.
In order to perform the integration, we must first “flip” the function p(t) to obtain p(−t), and then “shift” to obtain p(s − t). The discrete version of convolution is a summation over a finite number of terms. That is, pixels or values of the blurred image are obtained from a weighted sum of the corresponding pixel and its neighbors in the true image. The weights are given by the elements in the PSF. Consider a small example. Suppose discretizations of a true signal and PSF are given by
b(s) =
p(s − t)x(t)dt, 0 ≤ s ≤ 1 , (1)
w1 w2
x1 x2
x= x3 x4
x5 y 1
y2
p1
p2
and p = p3 (2)
p 4 p5
where wi and yi represent pixels in the original scene that are actually outside the field of view. The basic idea of computing the convolution, b, of x and p, can be summarized as follows:
• Rotate (i.e., flip) the PSF array p by 180 degrees, writing its elements from bottom to top.
• Match coefficients of the rotated PSF array with those in x by placing (i.e., shifting) the center of the PSF array (i.e., the entry corresponding to a shift of zero) over the ith entry in x.
• Multiply corresponding components, and sum to get the ith entry in b.
1The book Deblurring Images: Matrices, Spectra, and Filtering by Hansen, Nagy and O’Leary, 2006 is a good reference for image deconvolution.
For example, assuming that p3 is the center of the PSF array, we have
b1 =p5w1 +p4w2 +p3x1 +p2x2 +p1x3 b2 =p5w2 +p4x1 +p3x2 +p2x3 +p1x4 b3 = p5x1 +p4x2 +p3x3 +p2x4 +p1x5 b4 = p5x2 +p4x3 +p3x4 +p2x5 +p1y1 b5 = p5x3 +p4x4 +p3x5 +p2y1 +p1y2
1. Use the above equations to fill in the matrix-vector representation of 1D convolution
w1 w2
x1 b2 x2
b = b3 = x3 (3) b4 x4
b5 x5 y 1
y2
It is important to keep in mind that the values wi and yi contribute to the observed pixels in the blurred signal, even though they are outside the field of view. Since we do not know these values, boundary conditions are used to relate them to the pixels xi within the field of view.
Notice that matrix A has special structure. Before we get into different boundary con- ditions, let’s investigate this structure. An n × n matrix T is called a Toeplitz matrix if its entries are constant on each diagonal, and so the coefficients are often written as tij = ti−j. For example, a 5 × 5 Toeplitz matrix is usually represented as
t0 t−1 t−2 t−3 t−4
t1 t0t−1t−2t−3 T=t2 t1t0t−1t−2.
t3 t2 t1 t0 t−1 t4 t3 t2 t1 t0
Matrices with this kind of structure arise in many engineering applications. Notice that the matrix is completely defined by the first row and first column of T. Matlab has a function to create a Toeplitz matrix given the first column and first row. For example, the commands:
5
b1
c = [0 1 2 3 4];
r = [0 -1 -2 -3 -4];
T = toeplitz(c, r);
Page 2 of 7
will build the Toeplitz matrix:
If p is the 3rd column of T, then
p4 T=p5
On the other hand, if p is the 4th column of T, then
p4 p3 p2 p1 0
0 −1 −2 −3 −4
1 0−1−2−3 T=2 1 0 −1 −2.
3 2 1 0−1 43210
In image restoration, our computational problems involve n × n banded Toeplitz matrices
that can be defined by a single column vector of vector, p, is given by:
length n. In particular, suppose a column
.
0
p1 p2 p 3 p4
p3
p2 p3 p4 p 5
p1 p2 p3 p 4 p5
0
0 p1. p 2 p3
p1
p2 p=p3 p 4 p5
0
0 0
p5 p4 p3 p2 T=0 p5 p4 p3 0 0 p 5 p 4 0 0 0 p5
p1 p2. p 3 p4
Thus, to define T uniquely, we need to know where to put p in T.
2. Fill in the following function to create a banded Toeplitz matrix as defined above.
10
function T = bandToeplitz( p, k )
% Construct a banded Toeplitz where vector p is in the kth column
% Input:
% p – vector containing the entries of the kth column of T
% k – integer indicating the column of T
% Output:
% T – banded Toeplitz matrix
Page 3 of 7
This banded Toeplitz matrix is the basis from which we construct blurring matrices in image restoration. T may need to be modified slightly, depending on the type of boundary condition we impose; that is, in general we will need to construct matrices of the form
A=T+E
where E is an “extra” amount added to T which enforces certain boundary conditions. The
particular structure of E depends on the boundary condition:
• For zero boundary conditions where wi = yi = 0, E = 0.
• For periodic boundary conditions, E is defined so that A = T + E is a circulant matrix. For example:
p3 p2 p1 0 0 0 0 0 p5 p4
p4 p3 p2 p1 0 0 0 0 0 p5 If T=p5 p4 p3 p2 p1, then E=0 0 0 0 0
and
0 p 5 p 4 p 3 p 2 0 0 p5 p4 p3
p4 p3 p2 p1 0
p 1 0 0 0 0 p2 p1 0 0 0
0 0 0 0 p5
0 0 0 0 0 If T=0 p5 p4 p3 p2, then E=p1 0 0 0 0
p5 p4 p3 p2 p1
0 0 p5 p4 p3 p2 p1 0 0 0 0 0 0 p5 p4 p3 p2 p1 0 0
• For reflexive boundary conditions, E is a Hankel matrix. It is probably easiest to describe the exact form of E through some examples:
and
p3 p2 p1
p4 p3 p2 If T=p5 p4 p3 0 p 5 p 4 0 0 p5
p4 p3 p2
p5 p4 p3 If T=0 p5 p4 0 0 p 5
0 0 0
00 p4 p5000
p1 p2 p 3 p4
p1 p2 p3 p 4 p5
0 p5 0000 p1, then E=0 0 0 0 0 p 2 0 0 0 0 p 1 p3 0 0 0 p1 p2
0 p5
p1 0
p2, then E=0 p 3 0
p4 0
0000 0 0 0 0 0 0 0 p1 0 0 p 1 p 2 0 p1 p2 p3
Page 4 of 7
20
3. Next fill in the following Matlab function that will create the matrix E that results from the imposed boundary condition.
function E = buildBC(T, BC)
% This function builds the extra bit that needs to be
% added to the Toeplitz matrix to compensate for the
% boundary condition in image restoration problems.
% Input:
switch BC
case ‘zero’
E = zeros(size(T));
case ‘periodic’
E = buildPeriodic(T);
case ‘reflexive’
E = buildReflexive(T);
otherwise
error(‘Boundary condition is not valid.’)
end
function E = buildPeriodic(T)
banded Toeplitz matrix
character string specifying the desired boundary condition.
Choices are: ‘zero’, ‘periodic’, and ‘reflexive’
%
%
%
% Ouput:
% E – matrix having the same dimension as T.
%
T- BC-
% % % %
This subfunction imposes a periodic boundary condition by
creating a matrix, E, the extra bit that needs to be
added to the Toeplitz matrix, T, to create A = T + E,
which is a circulant matrix.
function E = buildReflexive(T)
% % %
This subfunction imposes a reflexive boundary condition by
creating a Hankel matrix, E, the extra bit that needs to be
added to the Toeplitz matrix, T, to create matrix A = T + E.
That is, you need to write sub-functions buildPeriodic and buildReflexive. Notice that these sub-functions have the same form as a normal function, but in this case, can and should be placed in the same m-file as the function buildBC. In buildPeriodic, you should use Matlab functions fliplr, flipud, and toeplitz to construct E using at most 3 lines of code. For buildReflexive you only need the Matlab function hankel, and again, you should be able to construct E in at most 3 lines of code.
Now that you have written the codes to generate the blur matrix for various boundary conditions, we can generate a 1D convolution example. Download the file blur1D.m from Canvas. This code calls your functions bandToeplitz.m and buildBC.m to generate a signal blurring example.
Page 5 of 7
5
4. In the file blur1D.m, the user can select from three example signals provided. The first one represents a bar code signal. The second is a smooth signal and the third is an edgy signal. Find an application for signal deconvolution and create an additional test signal that is representative of desired signals in your application. Include your signal in blur1D.m as example case 4. Explain why you created this particular test signal (i.e., why did you form the signal in this way and what do the features represent?). Feel free to be creative in your choice of application.
5. The next problem investigates the input parameter σ to vary the amount of blur. For one signal, construct blurred signals for values of σ = 1,2,3 and provide all blurred signals on one plot, along with the true signal.
6. Here we investigate the effect of different boundary conditions on the convolution matrix A. For fixed σ, generate convolution matrices corresponding to zero, periodic, and reflexive boundary conditions. Use matlab function spy and subplot to display the structures of A in one figure.
7. The purpose of this exercise is to illustrate how the SVD can be used to analyze the smoothing effects of a first-kind Fredholm integral equation. We use the one-dimensional reconstruction test problem, which is implemented in Regularization Tools as function shaw. The kernel in this problem is given by
2 sin(π(sin(s) + sin(t)))2 K(s, t) = (cos(s) + cos(t)) π(sin(s) + sin(t)) ,
−π/2 ≤ s, t ≤ π/2, while the solution is
f (t) = 2 exp(−6(t − 0.8)2 ) + exp(−2(t + 0.5)2 ).
This integral equation models a situation where light passes through an infinitely long slit, and the function f(t) is the incoming light intensity as a function of the incidence angle t. The problem is discretized by means of the midpoint quadrature rule to produce A and xexact, after which the exact right-hand side is computed as bexact = Axexact. The elements of bexact represent the outgoing light intensity on the other side of the slit.
(a) Choose n = 24 and generate the problem. Then compute the SVD of A and provide a plot of the left and right singular vectors. You should use the subplot command to create a 4 × 6 grid of plots, where each plot contains the vectors ui and vi. Describe your observations about the number of sign changes in these vectors as i gets larger?
(b) Use the function picard from Regularization Tools to inspect the singular values σi and the absolute SVD coefficients |ui⊤bexact|, as well as the corresponding absolute solution coefficients |ui⊤bexact|/σi. Provide the Picard plot. Is the discrete Picard condition satisfied? Explain.
10
10
40
Page 6 of 7
(c) Add a very small amount of noise e to the right-hand side bexact, b = bexact + e with ∥e∥2 / ∥bexact∥2 = 10−10. Provide the Picard plot again. What happens to the absolute SVD coefficients |ui⊤b| corresponding to the small singular values?
(d) Recall that the undesired “naive” solution x = A−1b can be written in terms of
the SVD as
Compute the partial sums
x = A−1b = n ui⊤bvi. i=1 σi
xk=k ui⊤bvi i=1 σi
for k = 1, 2, …. Use the subplot command to create a 4 × 6 grid of plots, where the each plot contains the vector xk along with the true solution for k = 1, 2, …, 24. Explain the behavior of these vectors.
Be sure to label all figures. You should become very familiar with these codes, since we will be using them as test examples in later classes.
Page 7 of 7