联系方式

  • QQ:99515681
  • 邮箱:99515681@qq.com
  • 工作时间:8:00-23:00
  • 微信:codinghelp

您当前位置:首页 >> Matlab编程Matlab编程

日期:2019-12-14 10:51

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


版权所有:留学生编程辅导网 2020 All Rights Reserved 联系方式:QQ:99515681 微信:codinghelp 电子信箱:99515681@qq.com
免责声明:本站部分内容从网络整理而来,只供参考!如有版权问题可联系本站删除。 站长地图

python代写
微信客服:codinghelp