CSED551 –HOMEWORK 4 solved

Description

.
NON-BLIND DECONVOLUTION
Image deblurring is a long-standing image restoration problem in computer vision and image processing. Most
classical image deblurring methods model blur using a convolution operation. Specifically, they model a blurred image
𝑏 as:
𝑏 = 𝑘 ∗ 𝑙 + 𝑛
where 𝑘 is a blur kernel or a point spread function (PSF), 𝑙 is a latent sharp image, and 𝑛 is additive noise, which
is often assumed to be white Gaussian such that 𝑛~𝑁(0, 𝜎
2
) where `white’ noise means the noise variance is constant
for all frequency bands as will be explained later. * is the convolution operator. Then, image deblurring becomes a
deconvolution problem, in which we invert the convolution operation. If we know both 𝑏 and 𝑘, the problem is
called non-blind deconvolution. On the other hand, if we know only 𝑏, the problem is called blind deconvolution.
Non-blind deconvolution serves as a fundamental building block in many applications, e.g., removing motion blur,
out-of-focus blur, and chromatic aberration.
INVERSE FILTERING
The most naïve approach to non-blind deconvolution is the inverse filtering. We can derive the inverse filtering by
applying the Fourier transform to the blur model above. Specifically, we can obtain:
𝐵 = 𝐾 ⊗ 𝐿 + 𝑁
where 𝐵 is the Fourier transform of 𝑏, i.e., 𝐵 = ℱ(𝑏) where ℱ represents the forward Fourier transform function.
𝐾, 𝐿, and 𝑁 are also the Fourier transforms of 𝑘, 𝑙, and 𝑛, respectively. ⊗ is the element-wise multiplication
operator. Ignoring noise 𝑁, we can obtain 𝐿 as:
𝐿 ≈
𝐵
𝐾
where 𝐵
𝐾
is the pixel-wise division between 𝐵 and 𝐾 in the frequency domain. Then, by computing the inverse
Fourier transform, we can obtain a deconvolution result. As this technique directly applies the inverse of a blur kernel
𝑘 to the input blurred image 𝑏 in the frequency domain, it is called the inverse filtering.
The inverse filtering, however, cannot produce satisfying results. The blur kernel 𝐾 in the frequency domain may
have elements whose values are zero or close to zero, and dividing 𝐵 with such elements leads to totally corrupted
images with significantly amplified noise. To overcome this issue, we need proper handling of noise and regularization.
WIENER DECONVOLUTION
One classical approach to non-blind deconvolution is the wiener deconvolution. The wiener deconvolution is derived
by minimizing the expected error in the frequency domain. Specifically, we want to find a deconvolution filter 𝐻 in
the frequency domain such that
min
𝐻
𝐸[‖𝐿 −𝐻𝐵‖
2
]
where 𝐻𝐵 is the pixel-wise multiplication between 𝐻 and 𝐵 in the frequency domain. By solving the equation
above, we can obtain a deconvolution filter 𝐻 such that 𝐻𝐵 is a deconvolution result with the least expected squared
error.
By replacing 𝐵 in the equation above with 𝐵 = 𝐾𝐿 + 𝑁, we can obtain:
min
𝐻
𝐸[‖𝐿 − 𝐻(𝐾𝐿 +𝑁)‖] = min
𝐻
𝐸[‖(1 −𝐻𝐾)𝐿 − 𝐻𝑁‖]
Then, we can extend the squares and obtain:
min
𝐻
𝐸[‖(1 − 𝐻𝐾)𝐿 − 𝐻𝑁‖] = min
𝐻
{‖1 −𝐻𝐾‖
2𝐸[‖𝐿‖
2
] −2(1 − 𝐻𝐾)𝐸[𝐿𝑁] + ‖𝐻‖
2𝐸[‖𝑁‖
2
]}
To further continue the derivation, we assume that noise 𝑁 is zero-mean Gaussian noise, which is independent to 𝐿
(in practice, it’s not, but many image restoration works use this assumption to simplify their problems.). Then
𝐸[𝐿𝑁] = 𝐸[𝐿]𝐸[𝑁]. As 𝑁 is zero-mean, 𝐸[𝑁] = 0. As a result, we obtain:
min
𝐻
{‖1 − 𝐻𝐾‖
2𝐸[‖𝐿‖
2
] +‖𝐻‖
2𝐸[‖𝑁‖
2
]}
By differentiating the equation above with respect to 𝐻, and setting it to zero, we can solve the equation for 𝐻:
⇒ −2(1 − 𝐻𝐾)𝐸[‖𝐿‖
2
]+ 2𝐻𝐸[‖𝑁‖
2
] = 0
⇒ 𝐻 =
𝐾𝐸[‖𝐿‖
2
]
𝐾2𝐸[‖𝐿‖
2]+ 𝐸[‖𝑁‖
2]
=
𝐾
2
𝐾2 +
𝐸[‖𝑁‖
2]
𝐸[‖𝐿‖
2]
⋅
1
𝐾
Finally, multiplying 𝐻 to 𝐵, we can obtain the formulation of the wiener deconvolution:
𝑙 = ℱ
−1 (
|ℱ(𝑘)|
2
|ℱ(𝑘)|
2 +
1
𝑆𝑁𝑅(𝜔)
⋅
ℱ(𝑏)
ℱ(𝑘)
)
where all the multiplication and division operations are pixel-wise operations. 𝑆𝑁𝑅(𝜔) is the signal-to-noise ratio at
each frequency 𝜔, which is required to be estimated. 𝑆𝑁𝑅(𝜔) is defined as:
𝑆𝑁𝑅(𝜔) =
𝑠𝑖𝑔𝑛𝑎𝑙 𝑣𝑎𝑟𝑖𝑎𝑛𝑐𝑒 𝑎𝑡 𝜔
𝑛𝑜𝑖𝑠𝑒 𝑣𝑎𝑟𝑖𝑎𝑛𝑐𝑒 𝑎𝑡 𝜔
Natural images tend to have signal variance that scales as 1
𝜔2
. Noise variance tends to be a constant, i.e., 𝜎𝑁(𝜔) =
𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡 for all 𝜔. This is called white noise. As a result, we can assume
𝑆𝑁𝑅(𝜔) =
1
𝑐𝜔2
and obtain a simplified version of the wiener deconvolution as:
𝑙 = ℱ
−1 (
|ℱ(𝑘)|
2
|ℱ(𝑘)|
2 + 𝑐𝜔2
⋅
ℱ(𝑏)
ℱ(𝑘)
)
where 𝑐 is a constant. In practice, a more simplified version is often used, which is defined as:
𝑙 = ℱ
−1 (
|ℱ(𝑘)|
2
|ℱ(𝑘)|
2 + 𝑐
⋅
ℱ(𝑏)
ℱ(𝑘)
)
(1)
By setting 𝑐 = 0, the equation above reduces to the inverse filtering. In this homework, you need to implement this
simplified version of the wiener deconvolution. You also need to discuss the effect of the parameter 𝑐 with example
images.
TOTAL-VARIATION (TV) DECONVOLUTION
Non-blind deconvolution is an ill-posed problem. Due to the loss of information caused by the noise 𝑛, even if we
know 𝑘, there can be an infinite number of solutions for 𝑙. To resolve such ill-posedness, most methods adopt a
regularization term or a prior on the latent image 𝑙. One of the most widely-used prior is the total variation, which is
defined as:
𝐸TV(𝑙) = ‖∂x ∗ 𝑙‖1 +‖∂y ∗ 𝑙‖
1
where ‖𝑓‖1
is the 𝐿1 norm of an image 𝑓 which is defined as:
‖𝑓‖1 = ∑|𝑓(𝑥, 𝑦)|
𝑥,𝑦
and ∂x and ∂y are partial derivative filters along the x- and y-axes, e.g., ∂x = [0,−1,1] and ∂y = [0,−1,1]
T
. With
the prior 𝐸TV(𝑙), we can define an energy function for non-blind deconvolution as:
𝐸(𝑙) = ‖𝑏 − 𝑘 ∗ 𝑙‖
2 + 𝛼𝐸TV(𝑙)
where 𝛼 is a parameter to control the strength of the regularization term.
You can also formulate the non-blind deconvolution problem as a maximum-a-posteriori estimation problem, and
derive the energy function above from the MAP formulation. Specifically, you can define a posterior distribution:
𝑝(𝑙|𝑏, 𝑘) ∝ 𝑝(𝑏|𝑙, 𝑘)𝑝(𝑙)
where 𝑝(𝑏|𝑙, 𝑘) is a likelihood, and 𝑝(𝑙) is a prior distribution.
Optimization
The regularization term 𝐸TV(𝑙) in the energy function 𝐸(𝑙) is a non-quadratic term, so the optimization of 𝐸(𝑙)
needs a non-linear optimization algorithm such as the iterative reweighted least-squares (IRLS) method, which is
explained in the course material (see 15.global_optimization.pdf). At each iteration of the IRLS method, you need to
solve a least-squares problem, which is formulated as:
𝐥
∗ = arg min
𝐥
{‖𝐛 −𝐊𝐥‖
2 + 𝛼(𝐥
T𝐃𝑥
T𝐖𝑥𝐃𝑥𝐥+ 𝐥
T𝐃𝑦
T𝐖𝑦𝐃𝑦𝐥)}
where 𝐛 and 𝐥 are vector representations of 𝑏 and 𝑙, respectively. 𝐊 is a matrix representing the convolution
operation with 𝑘. The equation above can be re-written as:
𝐥
∗ = arg min
𝐥
{𝐥
T𝐊
T𝐊𝐥 −2𝐥
T𝐊
T𝐛 + 𝐛
T𝐛 + 𝛼(𝐥
T𝐃𝑥
T𝐖𝑥𝐃𝑥𝐥+ 𝐥
T𝐃𝑦
T𝐖𝑦𝐃𝑦𝐥)}
As the right-hand side of this equation is a quadratic function with respect to 𝐥, its solution can be found by setting its
derivative to zero. Specifically, the derivative of the right-hand side is derived as:
2𝐊
T𝐊𝐥− 2𝐊
T𝐛 +2𝛼(𝐃𝑥
T𝐖𝑥𝐃𝑥 + 𝐃𝑦
T𝐖𝑦𝐃𝑦)𝐥
By setting the derivative to zero, we can derive the following equation:
(𝐊
T𝐊 + 2𝛼(𝐃𝑥
T𝐖𝑥𝐃𝑥 + 𝐃𝑦
T𝐖𝑦𝐃𝑦)) 𝐥 = 𝐊
T𝐛
Note that this equation has the same form as 𝐀𝐱 = 𝐛, which means that you can find a solution 𝐥 using a linear solver.
However, the size of the matrix (𝐊
T𝐊+ 2𝛼(𝐃𝑥
T𝐖𝑥𝐃𝑥 +𝐃𝑦
T𝐖𝑦𝐃𝑦)) is huge, so it is hard to directly solve the linear
system. Instead, you can use an iterative solver such as the conjugate gradient method. For more details, refer to the
reference code that is provided with this homework document.
REQUIREMENTS
In this homework, you need to implement two non-blind deconvolution algorithms that remove blur from a given
blurred image 𝑏 and a blur kernel 𝑘. Specifically, you are required to:
⚫ Wiener deconvolution
◼ Implement the wiener deconvolution according to Equation (1).
◼ Discuss the effect of the parameter 𝑐 with example images.
⚫ TV deconvolution
◼ Formulate the non-blind deconvolution problem as a maximum-a-posteriori (MAP) problem, and
derive the energy function 𝐸(𝑙) described above.
◼ Explain the relationship between the MAP formulation, and the energy function. Specifically, you need
to explain how 𝛼 can be derived from your MAP formulation.
◼ You need to implement an algorithm to optimize the energy function. The energy function above
involves 𝐸TV(𝑙), which is non-quadratic. Thus, you need to implement a non-linear optimization
algorithm. In this homework, you can implement the iterative-reweighted-least-squares (IRLS) method,
which was presented in our course material. See the file ’15.global_optimization.pdf’.
The IRLS optimization for non-blind deconvolution needs to compute convolution operations with a blur kernel
several times, which can take an excessive amount of computation time for large blur kernels. To accelerate the
computation, you can use the convolution theorem.
Your report must include:
⚫ Derivation of the energy function from the MAP formulation and discussion about their relationship
⚫ Detailed discussion on your implementation of the wiener deconvolution and TV deconvolution
⚫ Your results with a detailed discussion including the discussion on the parameter 𝑐 of the wiener
deconvolution
⚫ Limitations of the techniques that you found
You must upload a single zip file that contains the following to the LMS:
⚫ code/ – a directory containing all your code for this assignment
⚫ images/ – a directory containing your input images and their results
⚫ report.pdf – your report as a PDF file
Due: Dec. 19
th, 23:59
Penalty for late submission
⚫ 1 day: 70%
⚫ 2 days: 30%
⚫ 3 days: 0%
RUBRIC
⚫ 30 pts: wiener deconvolution
⚫ 20 pts: Energy function derivation for the TV deconvolution
⚫ 30 pts: Implementation of the TV deconvolution
⚫ 20 pts: Report
Your homework will be scored based on your report, and I am not going to compile or run your code. Thus, your
report must include all necessary details of your implementation and results.