|
MDS6106作业代做、代做Decision Analytics作业、代写MATLAB实验作业、MATLAB编程语言作业调试
iDDA · Institute for Data and Decision Analytics
MDS6106 – Introduction to Optimization
Final Project
Image Inpainting and Nonconvex Image Compression
The goal of this project is to investigate different optimization models and to utilize minimization
methodologies that were introduced in the lecture to reconstruct images from partial data.
Project Description. Typically, a grey-scale image U ∈ R
m×n
is represented as a m × n matrix
where each entry Uij represents a specific pixel of the image containing the color information. If
the columns of U = (u[1], u[2], . . . , u[n]) are stacked, we obtain the vector form
u = vec(U) = (u>[1], u>[2], . . . , u>[n])> ∈ R
mn
of the image U. In this project, we consider a class of imaging problems that is known as inpainting
problems. In general, inpainting describes the task of recovering an image U from partial data
and observations. Specifically, we assume that parts of the image U are missing or damaged, (e.g.,
due to scratches, stains, or compression), and the aim is to reconstruct this missing or damaged
information via solving a suitable inpainting or optimization problem.
Figure 1: Examples of different damaged images. We want to recover the missing image information
in the red areas.
The images in Figure 1 show several typical situations where inpainting techniques can be applied.
Our overall task is to reconstruct the red target areas. In this project, we assume that these
target areas are known, i.e., we have access to a binary mask Ind ∈ R
m×n with
Indij =
(
1 the pixel (i, j) is not damaged,
0 the pixel (i, j) is damaged.
This mask contains the relevant information about missing or damaged parts in the image.
Let us set ind := vec(Ind) ∈ R
mn
, s := Pmn
i=1 indi
, and Ω := {i : indi = 1}. Moreover, let
Ω = {q1, q2, ..., qs} denote the different elements of the index set Ω. Then, we can define the
following selection matrix
The vectors ej are unit vectors in R
mn and the matrix P selects all undamaged pixels of a stacked
image according to the mask ind. In particular, if U ∈ R
m×n
is the input image and u = vec(U)
is its stacked version, then the vector
b = P u ∈ Rs
contains the color information of all undamaged pixels of the original image U. Hence, we now
want to find a reconstruction y ∈ R
mn of u such that:
(i) Original information of the undamaged parts in U is maintained, i.e., x satisfies
P y = b or P y ≈ b. (1)
(ii) The image y can recover the missing parts in U in a “suitable” way.
In this project, we will discuss different models and strategies to achieve these goals.
Project Tasks.
1. Sparse Reconstruction and L-BFGS. The linear system of equations (1) is underdetermined
and has infinitely many possible solutions. In order to recover solutions with a “natural
image structure”, we first consider the so-called `1-regularized image reconstruction problem
is derived as in (1) and µ > 0 is given. The efficiency and success of this
model is based on the fact that the `1-regularization promotes sparsity and that there exist
canonical sparse representations of the image data. The idea is to use a linear transformation
x = Ψy of the image y in (1) and to transfer it to the frequency domain. In this new space,
many images have a very sparse representation and many of the components xi are zero or
close to zero. This motivates the choice of the `1-norm, kxk1 =
Pmn
i=1 |xi
|, in the model (2).
The quadratic term in (2) corresponds to the condition (1) and is small when the pixels of
undamaged parts of u and of the corresponding reconstruction have close values.
In this part, we want to use the discrete cosine transformation as sparse basis for the images,
i.e., we can set x = Ψy = dct(y). Since the DCT-transformation is orthogonal, this leads to
the following definitions:
Ax = A(x) = Pidct(x), A>v = A
>(v) = dct(P
>v), x ∈ R
mn, v ∈ R
s,
where idct denotes the inverse DCT-transformation. In MATLAB, these operations can be
represented compactly as functions via
P = @(x) x(ind); A = @(x) P(idct(x)); A
>A = @(x) dct(ind.*idct(x)).
Since the `1-norm is not differentiable, we also want to consider a smoothed and popular
variant of the `1-problem: (3)
where the so-called Huber-function is defined as follows:
ϕhub(z) = ( 12δz2if |z| ≤ δ,|z| − 12δ if |z| > δ,z ∈ R, δ > 0.
In contrast to the `1-norm, the Huber function is continuously differentiable and hence,
gradient-type methods can be applied to solve the problem (3).
2 of 6
– Implement the accelerated proximal gradient method (FISTA) for the `1-problem (2).
(It holds that λmax(A>A) ≤ 1).
– Implement the globalized L-BFGS method (Algorithm 5 with L-BFGS updates) for
the smooth Huber-loss formulation (3). You can use backtracking to perform the line
search and the two-loop recursion to calculate the L-BFGS update. You can choose
H0
as initial matrix for the update. In order to guarantee positive definiteness of the
L-BFGS update, the pair {s
k
, yk} should only be added to the current curvature pairs if
the condition (sk)>yk > 10−14 is satisfied. Suitable choices for the memory parameter
are m ∈ {5, 7, 10, 12}.– A suitable image quality measure is the so-called PSNR value. Suppose that u∗ =vec(U∗) is the original true (and undamaged) image, then we have:
PSNR := 10 · log10 mn
ky − u∗k2
where y = idct(x).
Compare the results of the two methods and models with respect to the PSNR value,
the number of iterations, and the required cpu-time. Test different parameter settings
and variants of the methods to improve the performance of your implementations. Plot
the reconstructed images and the convergence behavior of FISTA and L-BFGS. You
can use the stopping criterion kx
k+1 − y
k+1k ≤ tol in FISTA. How does the PSNR
value behave for different tolerances?
– The choice of the parameter µ and δ can depend on the tested images; good choices
are µ ∈ [0.001, 0.1] and δ ∈ [0.01, 0.5]. (The Huber function is closer to the absolute
value | · | for smaller choices of δ).
Hints and Guidelines:
– On Blackboard, we have provided two datasets containing test images and test inpainting
masks. Compare the performance of your methods on several images and different
type of masks.
– Typically, the pixel values Uij are represented by integer numbers in [0, 255]. Scale the
images to [0, 1]. A common initial point is zero: x
0 = zeros(m*n,1).
– For debugging, you can test your implementation of the L-BFGS strategy first on a
simpler example.
2. Total Variation Minimization. The total variation minimization problem utilizes a different
regularization to solve the inpainting task (1). The optimization problem is given by
min
x
ϕ(Dx) s.t. kP x − bk2 ≤ δ, (4)
where P ∈ R
s×mn and b ∈ R
s are given as in (1) and δ > 0 is a noise parameter. Here,
D ∈ R
2mn×mn is a discretized version of the image gradient and the mapping ϕ : R
2mn → R
is given by
ϕ(z) := Xmn
For an input image x = vec(X) and at the pixel (i, j), the image gradient Dx computes the
difference between the pixel values
δ1 = X(i+1)j − Xij and δ2 = Xi(j+1) − Xij
3 of 6
in the x- and y-direction of the image. The total variation term ϕ(Dx) penalizes now the
`2-norm of (δ1, δ2)
> at all pixels. In particular, the objective function in (4) is minimized,
when neighboring pixels have similar values.
In contrast to the `1-problem, this regularization also works in the pixel domain and we do
not need to transform the image x to the frequency space.
– Solve the optimization problem (4) and implement the primal-dual method. A pseudocode
of the primal-dual method can be found below (the algorithm will be discussed in
detail in the following lectures). You can use the estimate: kDk
2 ≤ 8.
– Run your code for different test images and masks and report your results. Repeat
your experiments with different step sizes σ, τ > 0 – can different choices improve the
performance of your implementation? Compare the PSNR values with the results in
part 1.) – does the total variation model achieve better results? The parameter δ can
be chosen small: δ ∈ {10−6
, 10−4
, 10−2}.
Hints and Further Guidelines:
– On BB, we have provided a MATLAB function D = get_D(m,n) to compute the operator
D for given dimensions m and n. You can use this code (or your own solutions) to
generate D. Note that the code is tailored to the definition of the function ϕ.
– The primal-dual method will be discussed in Lecture 23, November 27th in detail. It is
designed for problems of the form
min
x
f(x) + g(Kx),
where f : R
n → R and g : R
m → R are convex functions or indicator functions of
convex, closed sets, and K ∈ R
m×n
is a given matrix. The method is defined as follows:
• Initialization: Choose initial points x
0 ∈ Rn, y0 ∈ Rm and set x¯0 = x0. Choose
step sizes τ, σ > 0.• For k = 0, 1, . . . : Repeat the following steps:
Convergence of this method can only be guaranteed if the step sizes τ and σ are chosen
such that τσkKk
2 < 1. The primal-dual method does not have a simple or natural
stopping criterion. In practice, the algorithm is terminated after a suitable number of
iterations or when the relative change in xk
is small:
3. Nonconvex Image Compression. In this final part of the project, we try to investigate a
slightly different question related to inpainting. Given an image u ∈ R
mn, can we construct
a binary mask c = ind ∈ {0, 1}
mn such that s =
Pmn
i=1 ci
is small as possible, but the image
u can still be reconstructed with reasonable quality by only using pixels j ∈ {1, ..., mn} with
cj = 1? In this part, we want to consider a PDE-based image compression technique that
allows to compute such a mask c and the corresponding reconstruction x of u simultaneously.
The model is given by
min x,c12kx − uk2 + µkck1 s.t. diag(c)(x − u) − (I − diag(c))Lx = 0, c ∈ [0, 1], (5)
4 of 6
where u = vec(U) ∈ R
mn is the ground truth image, x ∈ R
mn is the reconstructed image,
c ∈ R
mn denotes the inpainting mask, and L = −D>D ∈ R
mn×mn is the Laplacian operator.
As long as one element in c is nonzero, the matrix A(c) = diag(c) + (diag(c) − I)L can be
shown to be invertible and hence, x can be expressed explicitly via x = A(c)
−1diag(c)u. In
this case, the problem (5) can be reduced to, this model has a standard form. Furthermore, the
gradient of f is given by
∇f(c) = diag(Lx + u − x)[A(c)>]−1(x − u) where x = A(c)−1diag(c)u.
In Figure 2, an exemplary output or solution of this problem is shown. The figure in the
middle depicts the mask c and the associated reconstruction x of u is shown on the right.
In this example, the mask has a density of s
mn ≈ 6.38% pixels, i.e., only 6.38% of the pixels
of the ground truth image (which is shown on the left) are used.
Figure 2: Image compression via optimal selection masks. Left: ground truth image u. Middle:
inpainting mask c. Right: reconstructed image x. The mask c only selects 6.38% of the
available pixels. The PSNR value of the recovered image is 35.1.
– Implement the inertial proximal gradient method (iPiano) that was presented in the
lecture to solve the nonconvex optimization problem (6).
You can use the method get_D to generate the Laplacian operator L = −D>D. We
suggest to use c
0 = ones(m*n,1) as initial point and βk ≡ β = 0.8 (or a different value
close to 0.8). As the Lipschitz constant ` of the gradient ∇f is unknown, the simple
line search-type procedure of the inertial proximal gradient method can be used to
determine ` adaptively.
Hints and Guidelines:
– The parameter µ depends on the specific ground truth image u. Possible choices are
µ ∈ [0.001, 0.01]. If µ is too large, the mask c will be set to zero which will cause
numerical issues and the reconstruction fails.
– This problem can be computationally demanding. Try to implement iPiano as efficiently
as possible! You can test your implementation first on smaller images. For instance
the MATLAB-command imresize(u,0.5) can be used to scale the original image u to
half of its size.
5 of 6
– Ensure that the matrix A(c) is in a sparse format. In this case, the backslash operator
or similar methods for sparse linear systems can be utilized to compute the solution
x = A(c)
−1diag(c)u efficiently.
– In order to prevent the Lipschitz constant ` from being too large, you can decrease it
by a certain factor after a fixed number of iterations, i.e., a possible strategy is
if mod(iter,5) == 0, ` = 0.95 · `, α = 1.99(1 − β)/`, end
You can also experiment with continuation strategies for the parameter µ: start with
a larger parameter µ0 > µ and reduce it after a fixed number of iterations until it
coincides with the original parameter µ.
– You can either use a number of iterations (250, 500, or 1000, . . .) or the condition,
as stopping criterion.
Project Report and Presentation. This project is designed for groups of four students. Please
send the following information to 217012017@link.cuhk.edu.cn or andremilzarek@cuhk.edu.cn
until November, 29th, 6:00 pm:
• Name and student ID number of the participating students in your group, group name.
Please contact the instructor in case your group is smaller to allow adjustments of the project
outline and requirements.
A report should be written to summarize the project and to collect and present your different
results. The report itself should take no more than 15–20 typed pages plus a possible additional
appendix. It should cover the following topics and items:
• What is the project about?
• What have you done in the project? Which algorithmic components have you chosen to
implement? What are your main contributions?
• Summarize your main results and observations.
• Describe your conclusions about the different problems and methods that you have studied.
You can organize your report following the outlined structure in this project description. As the
different parts in this project only depend very loosely on each other, you can choose to distribute
the different tasks and parts among the group members. Please clearly indicate the responsibilities
and contributions of each student and mention if you have received help from other groups, the
teaching assistant, or the instructor.
Try to be brief, precise and fact-based. Main results should be presented by highly condensed and
organized data and not just piles of raw data. To support your summary and conclusions, you
can put more selected and organized data into an appendix which is not counted in the page limit.
Please use a cover page that shows the names and student ID numbers of your group members.
The deadline for the report submission is December, 16th, 15:00 pm. Please send your report
and supporting materials to 217012017@link.cuhk.edu.cn.
The individual presentations of the project are scheduled for December, 17th, afternoon or
December, 18th. More information will follow here soon.
6 of 6
因为专业，所以值得信赖。如有需要，请加QQ：99515681或邮箱：99515681@qq.com 微信：codehelp
|
