微信客服：xiaoxionga100
微信客服：ITCS521
r代写-MATH3609
时间：2021-04-11
MATH3609 Optimization, Networks and Graphs:
Coursework 1 and Coursework 2
Dr Nathan Broomhead and Dr Julian Stander
Module Taught in Semester 2, Academic Year 2020/21
1 General Information about Coursework 1 and
Coursework 2
Please read the following points before attempting the courseworks:
• The deadline for both Coursework 1 and Coursework 2 is 10 am on Thursday, 22nd
April, 2022. You should submit Coursework 1 and Coursework 2 separately through
the MATH3609 Optimization, Networks and Graphs DLE site. Your submission will be
marked anonymously.
• Coursework 1 and Coursework 2 are Group Courseworks. Please work in
a self-assigned group of up to four people. You will need to register your
group on the DLE in the usual way.
Within each coursework, each member of the group will receive the same mark, unless
any member chooses to make use of the Peer Assessment option. You should keep notes
of all your group meetings to use as evidence in case you choose to make use of the Peer
Assessment option. If you wish to make use of the Peer Assessment option, you will need
to make an appointment to see Dr Nathan Broomhead or Dr Julian Stander.
• You should give due consideration to your personal time management to ensure that
your work is submitted in plenty of time prior to the deadline. For further information,
please consult the Student Handbook
https://www.plymouth.ac.uk/your-university/governance/student-handbook,
https://www.plymouth.ac.uk/student-life/your-studies/essential-information/exams/exam-rules-and-
regulations/extenuating-circumstances,
https://www.plymouth.ac.uk/uploads/production/document/path/18/18487/Extenuating_Circumstances_
Policy_and_Procedures_20-21_v.3.pdf and
https://liveplymouthac.sharepoint.com/sites/x70/SitePages/Student-QA.aspx#q-a-session-1-safety-net-
video%E2%80%8B%E2%80%8B%E2%80%8B%E2%80%8B%E2%80%8B%E2%80%8B%E2%80%8B.
• Each member of the group can expect to spend up to around 15 hours on Coursework 1
and up to around 15 hours on Coursework 2. These suggested times do not include the
preparation work that you will need to do.
• Each coursework counts for 15% = 12 × 30% of your final mark on this module. Marks
will be assigned according to the marking grids on pages 26 and 30.
1
• Marked scripts will be returned within 20 working days of the submission date, and,
in particular, you will get full feedback on your work by Friday, 21st May, 2021, unless
extensions are used.
• You are reminded of the University’s Academic Regulations:
Academic offences occur when activity is undertaken which could confer an unfair
advantage to any candidate(s) in assessment. The University recognises the following
(including any attempt to carry out the actions described) as academic offences, regardless
of intent:
a. Plagiarism, which is copying or paraphrasing of other people’s work or ideas into a
submitted assessment without full acknowledgement. More information on plagiarism
is available here:
https://www.plymouth.ac.uk/student-life/your-studies/essential-information/regulations/plagiarism
b. Collusion, which is unauthorised collaboration of students (or others) in producing a
submitted assessment. The offence of collusion occurs if a student copies any part of
another student’s work, or allows their own work to be copied. Collusion also occurs
if other people contribute significantly to work that a student submits as their own.
The complete list of regulations can be found here:
https://www.plymouth.ac.uk/student-life/your-studies/essential-information/regulations
All group members automatically confirm that they have understood the University’s
policy on plagiarism and collusion when they submit each coursework.
We now state the relevant MATH3609 Optimization, Networks and Graphs Assessed Learning
Outcomes (ALOs) for this assignment.
At the end of the module the learner will be expected to be able to:
ALO1 Develop and analyse methods for the solution of optimisation problems;
ALO3 Demonstrate an understanding of the underlying concepts of graph theory;
ALO4 Use computer software for the solution of optimisation problems and report the
results.
You should keep these ALOs in mind when working on these courseworks.
2
2 Image Restoration: Preparation Work
You should study this section carefully. There has also been a lecture on it.
2.1 The Basic Idea
Let us begin by reading in an image and displaying it.
# Set the working directory
# This is where the image R_image_binary.csv is stored
#
# Please note that your working directory will be different
# from ours (only one of which is needed)
#
## setwd("/Users/julianstander/Dropbox/MATH3609_CW_2020_2021")
#
setwd("/Users/nbroomhead/Dropbox/MATH3609_CW_2020_2021")
#
R_logo <- as.matrix(read.table("R_image_binary.csv",
sep = ",")) # The image is thought of as a matrix
#
image(R_logo,
col = c("black", "white"))
0.0 0.2 0.4 0.6 0.8 1.0
0.
0
0.
2
0.
4
0.
6
0.
8
1.
0
Please note that this image is a binary image. It comprises two ‘colours’, black and white. In
Coursework 1 and Coursework 2 we work exclusively with binary images for simplicity.
We will discuss the axes labels later.
3
Unfortunately, often images are received in a degraded or noisy form such as
0.0 0.2 0.4 0.6 0.8 1.0
0.
0
0.
2
0.
4
0.
6
0.
8
1.
0
This could be because they are observed and then transmitted by a satellite and so are
subject to atmospheric distortions. Similarly, medical images are often distorted in some way;
just think of an ultrasound image of a baby or even an x-ray or magnetic resonance image
when there has been some patient movement.
4
The aim of image restoration is to try to restore the first image
0.0 0.2 0.4 0.6 0.8 1.0
0.
0
0.
2
0.
4
0.
6
0.
8
1.
0
from the degraded image
0.0 0.2 0.4 0.6 0.8 1.0
0.
0
0.
2
0.
4
0.
6
0.
8
1.
0
5
2.2 Some Mathematical Notation
Images are represented as matrices. Here is a small example to fix ideas:
x <- matrix(c(0, 0, 0,
0, 0, 1,
1, 1, 1,
0, 0, 1),
nrow = 4,
byrow = TRUE)
#
x
## [,1] [,2] [,3]
## [1,] 0 0 0
## [2,] 0 0 1
## [3,] 1 1 1
## [4,] 0 0 1
#
# Dimension
#
dim(x)
## [1] 4 3
nrow(x)
## [1] 4
ncol(x)
## [1] 3
#
# Specific
#
x[1, 2]
## [1] 0
x[3, 1]
## [1] 1
So the image that is stored in the R object x comprises 4 rows and 3 columns. The value
taken by the image in row 1, column 2 is 0 corresponding to black, while the value taken
by the image in row 3, column 1 is 1 corresponding to white. We say that the above image
comprises 4× 3 = 12 pixels. We shall use i to index the pixels: here i = 1, . . . , n = 12, where n
is the total number of pixels. The value taken at pixel i is denoted xi, where xi can take
the values 0 for black or 1 for white: xi ∈ {0,1}. This enables us to write the whole image as
the vector x = (x1, . . . , xn).
6
We label the pixels column by column:
Columns
R
ow
s
i = 1 or ( 1 , 1 )
i = 2 or ( 2 , 1 )
i = 3 or ( 3 , 1 )
i = 4 or ( 4 , 1 )
i = 5 or ( 1 , 2 )
i = 6 or ( 2 , 2 )
i = 7 or ( 3 , 2 )
i = 8 or ( 4 , 2 )
i = 9 or ( 1 , 3 )
i = 10 or ( 2 , 3 )
i = 11 or ( 3 , 3 )
i = 12 or ( 4 , 3 )
1 2 3
1
2
3
4
So that
x
## [,1] [,2] [,3]
## [1,] 0 0 0
## [2,] 0 0 1
## [3,] 1 1 1
## [4,] 0 0 1
is the vector
as.vector(x)
## [1] 0 0 1 0 0 0 1 0 0 1 1 1
or x = (x1, x2, . . . , x12) = (0,0,1,0,0,0,1,0,0,1,1,1).
7
Please note that in
Columns
R
ow
s
i = 1 or ( 1 , 1 )
i = 2 or ( 2 , 1 )
i = 3 or ( 3 , 1 )
i = 4 or ( 4 , 1 )
i = 5 or ( 1 , 2 )
i = 6 or ( 2 , 2 )
i = 7 or ( 3 , 2 )
i = 8 or ( 4 , 2 )
i = 9 or ( 1 , 3 )
i = 10 or ( 2 , 3 )
i = 11 or ( 3 , 3 )
i = 12 or ( 4 , 3 )
1 2 3
1
2
3
4
we also show the row and column numbers in brackets. We will not usually use these, but
they are added here for information. A point to notice is that the column labels increase
from left to right, while the row labels increase from top to bottom which is different from
standard Cartesian co-ordinates.
8
Let’s consider how R draws an image using its image function
x
## [,1] [,2] [,3]
## [1,] 0 0 0
## [2,] 0 0 1
## [3,] 1 1 1
## [4,] 0 0 1
#
image(x,
col = c("black", "white"))
0.0 0.2 0.4 0.6 0.8 1.0
−
0.
2
0.
2
0.
6
1.
0
#
# Or with better labels
#
n_row <- nrow(x) # Number of rows in x
n_row
## [1] 4
#
n_col <- ncol(x) # Number of columns in x
n_col
## [1] 3
#
image(1:n_row, # Label the rows
1:n_col, # Label the columns
x,
col = c("black", "white"),
xlab = "Row number",
9
ylab = "Column number",
axes = FALSE)
#
axis(1, 1:n_row)
axis(2, 1:n_col)
Row number
Co
lu
m
n
nu
m
be
r
1 2 3 4
1
2
3
So R thinks of the horizontal axis as representing row number, and the vertical axis as
representing column number, just like standard Cartesian co-ordinates: (row i, column j) =(horizontal co-ordinate x,vertical co-ordinate y). This means that the first column of the
matrix x is plotted horizontally at the bottom, the middle column of x is plotted horizontally
in the middle, and the final column of x is plotted horizontally at the top.
These plotting details are not important for what will come. However, we give them here for
information.
10
2.3 A Little More Mathematical Notation
We have seen that in practice we observe a degraded version of a ‘true’ image such as:
0.0 0.2 0.4 0.6 0.8 1.0
0.
0
0.
2
0.
4
0.
6
0.
8
1.
0
We refer to this observed or degraded image using the vector y = (y1, . . . , yn), where n is the
number of pixels. Here yi ∈ {0,1}, i = 1, . . . , n.
11
2.4 Some Literature and Other Irrelevant Information
In scientific writing, we should always give credit where credit is due. Many people have
contributed to statistical image restoration. The following are but two references among a
great many:
Besag, J. E. (1986) On the statistical analysis of dirty pictures (with discussion). Journal of
the Royal Statistical Society B, 48, 259–302.
Geman, D. and Geman, S. (1984) Stochastic relaxation, Gibbs distributions and the Bayesian
restoration of images. IEEE Transactions on Pattern Analysis and Machine Intelligence, 6,
721–41.
Coursework 1 and Coursework 2 are based on the work of
Greig, D. M., Porteous, B. T. and Seheult, A. H. (1989) Maximum a posteriori estimation
for binary images. Journal of the Royal Statistical Society B, 51, 271–279.
Greig, Porteous and Seheult (1989) assume black is coded as 1 and white as 0. In Coursework
1 and Coursework 2, we use take 1 to represent white and 0 to represent black, as this is
more natural and usual.
The work of Greig, Porteous and Seheult (1989) was expanded on in these works:
Jubb, M. D. (1989) Image Reconstruction. PhD thesis. University of Bath.
Stander, J. (1992) Some Topics in Statistical Image Analysis. PhD thesis. University of Bath.
Stander, J. and Silverman, B. W. (1994) Temperature schedules for simulated annealing.
Statistics and Computing, 4, 21–32.
Franconi, L. and Jennison, C. (1997) Comparison of a genetric algorithm and simulated
annealing in an application to statistical image reconstruction. Statistics and Computing, 7,
193–207.
No knowledge of any of this literature is assumed in either Coursework 1 or
Coursework 2.
We end this section by remarking that images are not usually binary. Often each pixel can
take a grey level so that xi ∈ {0,1, . . . ,255} where 0 is black and 255 is white for example,
or a colour, which is more complicated to represent. We work with binary images to keep
life relatively simple, as this may be the first time that you have met image reconstruction
problems.
12
2.5 Statistical Models
Image restoration proceeds by specifying statistical models for possible true images x =(x1, . . . , xn) and for degraded images y = (y1, . . . , yn).
Please note that this is not a statistics module, so not all the details that we present here are
needed for you to be able to do Coursework 1 or Coursework 2. We include these details
to explain how the image restoration problem can be expressed as an optimization
problem. It is with that optimization problem that we shall work.
2.6 A Statistical Model for y given x
We will assume the following statistical model for yi given xi:
Pr(yi = 0 ∣xi = 0) = 1 − p (no error made, no swap)
Pr(yi = 1 ∣xi = 0) = p (error made, swap from 0 to 1)
Pr(yi = 0 ∣xi = 1) = p (error made, swap from 1 to 0)
Pr(yi = 1 ∣xi = 1) = 1 − p (no error made, no swap)
Essentially, we are assuming that the colour at pixel i is swapped from black to white, or
from white to black, with probability p ∈ [0,1]; otherwise, it is not changed. For simplicity,
we assume that p is known.
We assume that the error mechanism at each pixel is independent of the error mechanisms at
every other pixel. This enables us to write down Pr(y ∣x):
Pr(y ∣x) = n∏
i=1 Pr(yi ∣xi).
Let’s recall that the degraded image y is observed, while the ‘true’ image x is unknown and
needs to be estimated. Thought of as a function of the unknown image vector x, Pr(y ∣x)
is the likelihood function that you may have met in MATH2602 Statistical Inference and
Regression or MATH3613 Data Modelling. We shall write L(y ∣x) = Pr(y ∣x).
13
Now let’s write f(yi ∣xi) = Pr(yi ∣xi), so that f(1 ∣0) = Pr(yi = 1 ∣xi = 0) for example. We
have
L(y ∣x) = Pr(y ∣x)
= n∏
i=1 Pr(yi ∣xi)
= n∏
i=1 f(yi ∣1)xif(yi ∣0)1−xi ,
since f(yi ∣1)xif(yi ∣0)1−xi = f(yi ∣1) if xi = 1 and f(yi ∣0) if xi = 0.
From this it follows that
L(y ∣x) = n∏
i=1 f(yi ∣0) × {f(yi ∣1)f(yi ∣0)}
xi
. (1)
Later, we will use λi, which is defined as follows:
λi = log{f(yi ∣1)
f(yi ∣0)} . (2)
Please note that λi only depends on the observed yi and not on xi. Please also note that, if
yi = 0, then
λi = log{f(yi ∣1)
f(yi ∣0)} = log{f(0 ∣1)f(0 ∣0)} = log ( p1 − p) ,
while, if yi = 1, then
λi = log{f(yi ∣1)
f(yi ∣0)} = log{f(1 ∣1)f(1 ∣0)} = log (1 − pp ) .
14
2.7 A Statistical Model for x
To define a statistical model for x, we need to explain what it means for two pixels to be
neighbours. This can be done by the following picture
Columns
R
ow
s
i = 1 or ( 1 , 1 )
i = 2 or ( 2 , 1 )
i = 3 or ( 3 , 1 )
i = 4 or ( 4 , 1 )
i = 5 or ( 1 , 2 )
i = 6 or ( 2 , 2 )
i = 7 or ( 3 , 2 )
i = 8 or ( 4 , 2 )
i = 9 or ( 1 , 3 )
i = 10 or ( 2 , 3 )
i = 11 or ( 3 , 3 )
i = 12 or ( 4 , 3 )
1 2 3
1
2
3
4
The red and blue lines indicate neighbouring pixles; we use blue to indicate horizontal
neighbours and red to indicate vertical neighbours, although this distinction is not important.
As an example we see that the neighbours of pixel i = 6 are 5, 10, 7 and 2. It turns out
that it simpler to assume that the image is wrapped on a torus, so that, for example the
neighbours of pixel i = 1 are 4 (wrapped top to bottom), 5, 2 and 9 (wrapped left to right).
This assumption has very little practical effect for large images, but it does mean that every
pixel has four neighbours.
15
These neighbourhood relationships lead to a symmetric adjacency matrix A, which here
will be of dimension n × n = 12 × 12:
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## [1,] 0 1 0 1 1 0 0 0 1 0 0 0
## [2,] 1 0 1 0 0 1 0 0 0 1 0 0
## [3,] 0 1 0 1 0 0 1 0 0 0 1 0
## [4,] 1 0 1 0 0 0 0 1 0 0 0 1
## [5,] 1 0 0 0 0 1 0 1 1 0 0 0
## [6,] 0 1 0 0 1 0 1 0 0 1 0 0
## [7,] 0 0 1 0 0 1 0 1 0 0 1 0
## [8,] 0 0 0 1 1 0 1 0 0 0 0 1
## [9,] 1 0 0 0 1 0 0 0 0 1 0 1
## [10,] 0 1 0 0 0 1 0 0 1 0 1 0
## [11,] 0 0 1 0 0 0 1 0 0 1 0 1
## [12,] 0 0 0 1 0 0 0 1 1 0 1 0
Do you understand why?
16
We can now define a probability mass function on images x as
Pr(x) = exp (−β × number of discrepant neighbouring pairs in x)∑images x′ exp (−β × number of discrepant neighbouring pairs in x′) ,
in which β ≥ 0 is a constant that we will discuss later. The sum over all images x′ comprises
2n terms, where n is the number of pixels, and can be very hard to compute when n is large.
As this sum does not depend on the argument x of Pr(x), we can discard it, if we replace
the equality with a ‘proportional to’ sign. So we write
Pr(x)∝ exp (−β × number of discrepant neighbouring pairs in x) .
For the image x that we have been working with:
## [,1] [,2] [,3]
## [1,] 0 0 0
## [2,] 0 0 1
## [3,] 1 1 1
## [4,] 0 0 1
the number of discrepant neighbouring pairs is 10:
Columns
R
ow
s
0
0
1
0
0
0
1
0
0
1
1
1
1 2 3
1
2
3
4
Hence, the probability of this image is Pr(x) = C exp(−β × 10), where C is a constant.
17
Similarly, the number of discrepant neighbouring pixels in the all black and the all white
images is, unsurprisingly, 0:
Columns
R
ow
s
0
0
0
0
0
0
0
0
0
0
0
0
1 2 3
1
2
3
4
Columns
R
ow
s
1
1
1
1
1
1
1
1
1
1
1
1
1 2 3
1
2
3
4
This means that the all black and the all white images are the most probable with probabilities
C exp(β × 0) = C each.
18
The chess-board image is the least probable, as it has the highest number of discrepant
neighbouring pixels:
Columns
R
ow
s
1
0
1
0
0
1
0
1
1
0
1
0
1 2 3
1
2
3
4
There are 20 discrepant pairs, meaning that the probability of the chess board image is
C exp(−β × 20).
You may be asking: what does β ≥ 0 do? The larger the value of β, the higher the relative
probabity given to the all black and to the all white image. In other words, higher β values
encourage ‘smoother’ images. Because of this, sometimes β is referred to as a smoothing
parameter.
Finally, we remark that our probability model for images x,
Pr(x)∝ exp (−β × number of discrepant neighbouring pairs in x) ,
is essentially the same as the Ising model used a lot in physics and statistical mechanics
to represent substances such as magnets. To quote Federico Camia: ‘The Ising model,
introduced by Lenz in 1920 to describe ferromagnetism, is one of the most studied models
of statistical mechanics. Its two dimensional version has played a special role in rigorous
statistical mechanics since Peierls’ proof of a phase transition in 1936 and Onsager’s derivation
of the free energy in 1944.’
19
It can be shown that the
number of discrepant neighbouring pairs in x = 12 n∑i=1 n∑j=1aij (xi − xj)2 ,
in which x = (x1, . . . , xn) is an image and aij is the (i, j)-th element of the adjacence matrix
A such that
aij = {1 if pixels i and j are neighbours0 if pixels i and j are not neighbours.
A consequence of this result is that
Pr(x) ∝ exp (−β × number of discrepant neighbouring pairs in x)
= exp(−β × 12 n∑i=1 n∑j=1aij (xi − xj)2)
= exp(−β2 n∑i=1 n∑j=1aij (xi − xj)2) . (3)
Another irrelevant fact is that more general probability mass functions can be defined by
replacing (xi − xj)2 with φ (xi − xj) where the function φ satisfies certain properties.
20
2.8 A Statistical Model for x given y
It makes sense to consider Pr(x ∣ y), as we are given the observed image y. Pr(x ∣ y) can be
found using Bayes’ theorem as follows:
Pr(x ∣ y) ∝ L(y ∣x) ×Pr(x)
∝ n∏
i=1 f(yi ∣0) × {f(yi ∣1)f(yi ∣0)}
xi ×Pr(x) from (1)
∝ n∏
i=1 f(yi ∣0) × {f(yi ∣1)f(yi ∣0)}
xi × exp(−β2 n∑i=1 n∑j=1aij (xi − xj)2) from (3) (4)
In the language of Bayesian statistical inference, Pr(x) is the prior probability mass
function of all images x before observing the degraded image y. Pr(x) expresses our belief
about x before seeing the data in terms of probabilities. Pr(x ∣ y) is our belief about x after
seeing the data y and is referred to as our posterior probability mass function.
Please note that, as mentioned, it does not matter if you do not follow all these details,
although they are not hard. The important part is that image restoration leads to an
optimization problem, as we shall see.
21
2.9 The Restoration of a Binary Image as an Optimization Problem
The restored image is defined as the most probable image under the posterior Pr(x ∣ y).
That is, it is the image x that maximizes Pr(x ∣ y). For this reason it is referred to as the
maximum a posteriori or MAP image.
From (4) we see that the restored image maximizes
Pr(x ∣ y)∝ n∏
i=1 f(yi ∣0) × {f(yi ∣1)f(yi ∣0)}
xi × exp(−β2 n∑i=1 n∑j=1aij (xi − xj)2)
over all images x.
We know that it’s very often easier to work with the logarithm of a function, rather than
the function itself, when maximizing. This is because the logarithm of a function and the
function itself are maximized at the same point, and also because taking logarithms turns
products into sums, which are easier to deal with. So, the restored image maximizes
log {Pr(x ∣ y)} = additive constant D + n∑
i=1 log [f(yi ∣0) × {f(yi ∣1)f(yi ∣0)}
xi] − β2 n∑i=1 n∑j=1aij (xi − xj)2
= D + n∑
i=1 log {f(yi ∣0)} + n∑i=1 log [{f(yi ∣1)f(yi ∣0)}
xi] − β2 n∑i=1 n∑j=1aij (xi − xj)2
= D + n∑
i=1 log {f(yi ∣0)} + n∑i=1 xi log{f(yi ∣1)f(yi ∣0)} − β2 n∑i=1 n∑j=1aij (xi − xj)2
= D + n∑
i=1 log {f(yi ∣0)} + n∑i=1 xiλi − β2 n∑i=1 n∑j=1aij (xi − xj)2 (5)
by the definition of λi in (2).
Please note that ∑ni=1 log {f(yi ∣0)} does not depend on x, and so can be considered a constant.
It therefore follows that we need to maximize
n∑
i=1 xiλi − β2 n∑i=1 n∑j=1aij (xi − xj)2
over x = (x1, . . . , xn), or, equivalently, to minimize
− n∑
i=1 xiλi + β2 n∑i=1 n∑j=1aij (xi − xj)2
over all images x.
22
So, we have finally arrived at our image restoration optimization problem. Given a degraded
image y = (y1, . . . , yn), we can find our restored image by solving the problem:
minimize: − n∑
i=1 xiλi + β2 n∑i=1 n∑j=1aij (xi − xj)2 over images x = (x1, . . . , xn). (6)
It can be shown that the minimiation problem given by (6) is equivalent to the problem
minimize:
n∑
i=1(1 − xi)max(0, λi) + n∑i=1 ximax(0,−λi) + β2 n∑i=1 n∑j=1aij (xi − xj)2 (7)
over images x, in which max(a, b) means the maximum of a and b, so that, for example,
max(0,1) = 1 and max(0,−1) = 0.
This quantity is never negative.
23
2.10 Solving the Minimization Problem using the Ford-Fulkerson
Algorithm
We will now see how Greig, Porteous and Seheult (1989) solved the minimization problem
(7), namely,
minimize:
n∑
i=1(1 − xi)max(0, λi) + n∑i=1 ximax(0,−λi) + β2 n∑i=1 n∑j=1aij (xi − xj)2
over images x, using the Ford-Fulkerson algorithm. To do this, Greig, Porteous and Seheult
(1989) defined a network with capacities as follows:
2.10.1 Definition of the Network with Capacities
• There are n + 2 vertices:
– A source s;
– A sink t;
– A vertex for each of the n pixels, indexed by i or j below.
• If λi > 0, there is an arc (directed edge) (s, i) from the source s to i with capactiy
csi = λi; otherwise
• If λi < 0, there is an arc (i, t) from i to the sink t with capactiy cit = −λi;
• There is an undirected edge {i, j} between two internal vertices i and j (that corre-
spond to pixels i and j) with capacities cij = β aij. This can be thought of as two arcs,
one between i and j, and the other between j and i, both of which have capacities β aij .
• There are no other arcs.
24
2.10.2 Definition of a Cut
Let us say that x = (x1, . . . , xn) is an image that is being proposed as a reconstruction of the
observed y = (y1, . . . , yn).
• Partitioning the Network
We can use x to partition the network into two sets:
W (x) = {s} ∪ {i ∶ xi = 1}, source s and all white pixels
B(x) = {i ∶ xi = 0} ∪ {t}, all black pixels and the sink t.
This partition of the vertices in the network into the sets W (x) and B(x) defines a cut.
• The Cut Capacity
The associated cut capacity is
C(x) ∶= C(∂+(W (x))) = ∑
kl∈∂+(W (x)) ckl. (8)
Notice that the arcs of ∂+(W (x)) are of three types: arcs starting at s; arcs ending at t; arcs
between intermediate vertices.
The cut capacity (8) can also be written as
C(x) = ∑
k∈W (x) ∑l∈B(x) ckl.
It turns out that
C(x) = n∑
i=1(1 − xi)max(0, λi) + n∑i=1 ximax(0,−λi) + β2 n∑i=1 n∑j=1aij (xi − xj)2 ,
which is exactly the function that we wanted to minimize in (7) to restore the image from y.,
So, the image restoration problem becomes a max-flow/min-cut problem: we find the cut
that minimizes the cut capacity C(x) and this determines x. We can do this using the
Ford-Fulkerson algorithm (or modern improvements of it). ,
We will expand on and explore the above in Coursework 1 and Coursework 2.
25
3 Coursework 1
Marks will be assigned for Coursework 1 according to this marking grid.
MATH3609 Optimization, Networks and Graphs:
Coursework 1 Marking Grid
Assessment Area Maximum Awarded Feedback
Mark Mark
Task 1: Correctness and va-
lidity of argument. Read-
ability of explanation.
15
Task 2: Correctness and va-
lidity of argument. Read-
ability of explanation.
20
Task 3: Correctness and va-
lidity of argument. Read-
ability of explanation.
10
Task4: Correctness and va-
lidity of argument. Read-
ability of explanation.
25
Task 5: Correctness and va-
lidity of argument. Read-
ability of explanation.
15
Task 6: Correctness and va-
lidity of argument. Read-
ability of explanation.
15
Total1 100
1This mark is provisional. Like all marks it is considered by an External Examiner and an Assessment Panel.
26
3.1 Coursework 1 Tasks
Coursework 1 comprises a number of tasks. You need to produce a document containing
your solutions to these tasks. These may be hand-written and then scanned or typed up
using, for example, LATEX. It is your responsibility to make sure that they are readable.
The page limit is fifteen pages. Please note that this is the limit and should not be
considered as the target. Well written, clear and concise mathematics will receive high credit.
Please arrange your report under subheadings referring to the tasks, that is, Task 1,
Task 2, . . .. This will help you to respond systematically to all parts of Coursework 1.
The marks for each task are shown in the grid on page 26.
Task 1: Let G = (V,E) be a graph with vertex set V = {1,2, . . . , n} and let A be the
adjacency matrix of G with (i, j)-th entry aij. Suppose x∶V → {0,1} is a function which
assigns to each vertex i a number xi, which is either 0 or 1. A pair of neighbouring vertices i
and j are called discrepant if xi ≠ xj.
Explain why the number of discrepant neighbouring (unordered) pairs is given by the formula:
1
2
n∑
i=1
n∑
j=1aij (xi − xj)2 .
Task 2: Explain why the minimization problem given by (6) on page 23 is equivalent to the
problem
minimize:
n∑
i=1(1 − xi)max(0, λi) + n∑i=1 ximax(0,−λi) + β2 n∑i=1 n∑j=1aij (xi − xj)2
over images x, in which max(a, b) means the maximum of a and b, so that, for example,
max(0,1) = 1 and max(0,−1) = 0.
Task 3: Prove that
n∑
i=1(1 − xi)max(0, λi) + n∑i=1 ximax(0,−λi) + β2 n∑i=1 n∑j=1aij (xi − xj)2
is never negative.
Task 4: Prove that
C(∂+(W (x))) = n∑
i=1(1 − xi)max(0, λi) + n∑i=1 ximax(0,−λi) + β2 n∑i=1 n∑j=1aij (xi − xj)2 .
Hint: The arcs of ∂+(W (x)) can be split into three types, namely, arcs that start at s, arcs
that end at t and arcs between intermediate vertices. Consider the sums of the capacities
arcs of each type separately.
27
Task 5: Show that every cut K of the network can be written as
K = ∂+(W (x))
for some image x, and deduce using this and the previous tasks that the Ford-Fulkerson
algorithm can be used to find a restored image.
Task 6: Is there always a unique solution x to the minimization problem (i.e. a unique
restored image)? Illustrate your answer with an example.
28
3.2 Coursework 1: What You Need to Submit
Your group needs to submit the following file electronically using the DLE:
• A Portable Document Format file MATH3609_Coursework_1_2020_2021.pdf containing
your report for Coursework 1.
All members of the group should take responsibility for seeing that your work
is submitted before the deadline.
It is also the responsibility of all members to check that the correct file has
been uploaded and that it is readable.
If anything is unclear, you should ask without delay.
29
4 Coursework 2
Marks will be assigned for Coursework 2 according to this marking grid.
MATH3609 Optimization, Networks and Graphs:
Coursework 2 Marking Grid
Assessment Area Maximum Awarded Feedback
Mark Mark
Image Restoration
Tasks: Defining the set-
up. Correctness, clarity
and presentation of code.
25
Image Restoration
Tasks: Optimizing,
presenting the results,
choosing β. Depth of
understanding and quality
of description of insights
gained.
50
Report: display results
clearly, concisely, system-
atically and professionally,
using correct spelling and
grammar. This is important
for employability.
25
Total2 100
2This mark is provisional. Like all marks it is considered by an External Examiner and an Assessment Panel.
30
4.1 Coursework 2 Tasks
Coursework 2 comprises a number of tasks. You need to produce a report of your work
following the instructions below.
Your report should contain well presented and annotated R code.
The page limit is twenty pages, including R code and figures. Please see page 39
for more details. Please do not submit an additional appendix as it will not be considered.
Reports that contain irrelevant or uninteresting discussion or R code will be penalized. Please
note that twenty pages is the limit and should not be considered as the target. Well written,
concise reports will receive high credit.
You may also omit lengthy, technical R code, if you feel that it does not add to
the argument that you are making. Essential code should, however, be shown.
Please arrange your report under subheadings referring to the tasks, that is, Task 1,
Task 2, . . .. This will help you to respond systematically to all parts of Coursework 2.
31
4.1.1 Implementing the Ford-Fulkerson Algorithm to Restore the
Image in R: Getting Started
Following on from the information in the section ‘Image Restoration: Preparation Work’ we
will now produce an image restoration from observed data y using the Ford-Fulkerson
algorithm implemented in R.
First, we need to read y into R. The data that we will use are supplied in the file
R_degraded_image_0.2.csv. We can read in and display these data in the same was
as we did at the beginning of this document. Here is the code again to make sure that you
get started properly:
y <- as.matrix(read.table("R_degraded_image_0.2.csv",
sep = ","))
#
image(y,
col = c("black", "white"))
0.0 0.2 0.4 0.6 0.8 1.0
0.
0
0.
2
0.
4
0.
6
0.
8
1.
0
#
n_row <- nrow(y)
n_col <- ncol(y)
#
c(n_row, n_col)
## [1] 23 28
#
y_vec <- as.vector(y)
#
length(y_vec)
## [1] 644
32
We can label the axes in a possibly nicer way:
image(1:n_row,
1:n_col,
y,
xlab = "Row number",
ylab = "Column number",
main = "Data: degraded image y",
col = c("black", "white"))
5 10 15 20
5
10
15
20
25
Data: degraded image y
Row number
Co
lu
m
n
nu
m
be
r
• It is not hard to imagine that y is a degraded version of the image in R_logo seen above.
• It is known that p = 0.2. (In practice, p and any other unknown parameter or parameters
can also be estimated form the data by assigning prior distributions to them and sampling
from the resulting posterior using Markov chain Monte Carlo methods as discussed in
MATH3613 Data Modelling. This is beyond the scope of Coursework 2.)
• We will assume for now that β = 1.
33
4.1.2 Defining the Vector of λ Values
Task 1:
Using the function ifelse or otherwise, define the vector of λ values (λ1, . . . , λn) in R and
show a few λ values. Identify and record in R which λi values are positive and which λi values
are negative using the which function, storing the results in vectors which_lambda_positive
and which_lambda_negative respectively.
4.1.3 Defining the Network with Capacities
We need to define a data frame comprising three variables or columns called from, to and
capacity. This data frame defines the capacities in the network. If, for example, we wish to
specify ckl, then the variable from would take the value k, the variable to would take the
value l and the variable capacity would take the value ckl.
One way to define the required data frame of capacities is to construct three data
frames df_source, df_sink and df_internal, and to combine df_source, df_sink and
df_internal by binding their rows into one data frame df using rbind.
The data frame df_source records the capacities of arcs coming from the source s.
Task 2: Define df_source in R. The square brackets [] may help you here. In order to
check this, the first few rows of the result are given
## from to capacity
## 1 source 2 1.386294
## 2 source 3 1.386294
## 3 source 4 1.386294
## 4 source 6 1.386294
## 5 source 7 1.386294
## 6 source 8 1.386294
The data frame df_sink records the capacities of arcs going to the sink t.
Task 3: Define df_sink in R and show a part of it.
34
Before we define df_internal, we need to obtain the adjacency matrix A. This can be
done by using the function get_neighbourhood_matrix with which you are supplied in the
file of helper functions helper_functions.R; you do not need to understand the details of
these helper functions. Please note that get_neighbourhood_matrix is defined using other
functions, so that you’ll have to make sure that all the functions in the file of helper functions
helper_functions.R are loaded into R. We also need to extract the rows and columns of A
corresonding to neighbouring pixels:
A <- get_neighbourhood_matrix(y, n_row)
#
A[1:5, 1:5]
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0 1 0 0 0
## [2,] 1 0 1 0 0
## [3,] 0 1 0 1 0
## [4,] 0 0 1 0 1
## [5,] 0 0 0 1 0
#
neighbouring_pixels <- which(A == 1,
arr.ind = TRUE)
#
head(neighbouring_pixels)
## row col
## [1,] 2 1
## [2,] 23 1
## [3,] 24 1
## [4,] 622 1
## [5,] 1 2
## [6,] 3 2
35
The data frame df_internal records cij where i and j correspond to pixels i and j.
Task 4: Define df_internal in R and show a part of it. Define the data frame df using
rbind. Use the function graph_from_data_frame from the igraph package to create an
igraph graph called g from df. Produce a plot of your capacitated network, as shown here;
please note, that it’s very big (but rather nice)! Please note that your plot may be displayed in
a different way from the one here. You may find https://igraph.org/r/doc/plot.common.html
helpful.
source
2
23
24
622
1
3
25
623
4
26
624
5
27
625
6
28
626
7
29
627
8
30
628
9
31
629
10
32
630
11
33
631
12
34
632
13
35633
14 36
634
15
37
635
16
38
636
17
39 637
18
40
638
19
41
639
20
42
640
21
43
64122
44
642
45
643
46
64447
48
49
50
51
52
53
54
55
56
57
58
5960
6162
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
8283
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124 125
126
127
128
129
130
131132
133
134
135
136
137
138
139
140
141
142
143
144
145
146 147
148
149
150
151
152153
154
155
156
157
158
159
160
161 162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195196
197198
199
200
201
202
203
204
205
206
207 208
209
210
211
212
213
214
215
216
217218219220
221
222
223
224
225
226
227
228
229
230 231
232 233
234
235
236
237
238
239
240
241242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264265
266
267
268
269
270
271
272
273
274 275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300 301
302
303
304
305
306
307
308
309
310
311
312313314
315
316
317
318
319
320
321 322
323
324
325
326
327
328
329330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384385
386
387
388
389
390
391
392
393
394
395
396
397
398
399400
401
402
40
404
405
406
407
408409
410
411
412
413
414
415
416
417
418
419
420
421
422423
424
425
426
427
428
429
430
431
432
433
434
435
436 437
438
439
440
441
442
443
444
445446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
47
472
473
474
475
476
477
478
79
480
481
482
483
484
485486487
488
489
490
491 492
493
494
495
496
497498
499
500
501
502
503
504
505
506
507
508
509
510
511
512513
514
515
516
517
518
519
520521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542543
544
545
546
547
548
549550
551
552553
554
555
556
557
558
559
560
561
562
563
564
565566567
568
569
570
571
572573574
575
576
577
578
579
580
581
582
583
584
585
586 587
588589590
591
592
593
594
595
596
597
598
599
600601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
sink
36
4.1.4 Producing the Image Restoration by Maximizing Flow or Minimizing Cut Capacity
The function max_flow provided by the igraph package finds the maximum flow in a graph
or equivalently the minimum cut capacity. max_flow effectively implements an enhanced
modern version of the Ford-Fulkerson algorithm. Amongst other output, max_flow
returns a partition of the vertices into two sets that it calls partition1 and partition2.
partition1 is the set W (x) and partition2 is the set B(x) defined above, corresponding to
the ‘maximum flow’/‘minimum cut’ x. partition1 and partition2 are quite complicated
R objects. However, the part of interest to us can be extracted using the names function;
a character vector is given by names. Once names has been applied and any unnecessary
element has been removed, the resulting character vector can be turned into pixels numbers
using the as.numeric function.
Task 5: Use max_flow to maximize the flow in g. Report the value of the maximum flow.
Working with partition1 find the numerical values of the white pixels in the set W (x),
corresponding to the max-flow/min-cut, and report the first six of these. Similarly, find the
numerical values of the black pixels in the set B(x) and report the first six.
Task 6: Use these pixel values to define the reconstructed image. An easy way (but not the
only way) to do this is to use rep to define a vector, x_map_vec say, of the same length as
y_vec comprising NAs, and then to assign the value 0 to pixels that have been identified as
black and the value 1 to pixels that have been identified as white. Then turn the vector into
a matrix, remembering that pixels are labelled column by column; please see the help file for
matrix, if necessary. Display your image. Here is the one that we obtained.
0.0 0.2 0.4 0.6 0.8 1.0
0.
0
0.
2
0.
4
0.
6
0.
8
1.
0
37
4.1.5 Choosing β
This reconstruction x obtained from the observed image y was obtained using β = 1. In fact,
the choice of β is in the hands of the user.
Task 7:
Discuss and illustrate how you could choose β.
There is no right answer here. What we are looking for is sensible suggestions for
the choice of β, discussed clearly and critically. Your discussion should be illus-
trated with useful figures, as appropriate. Concise answers comprising sensible
and well illustrated suggestions will receive high credit.
38
4.2 Coursework 2: Report Production
The ability to report technical results clearly and concisely is an important skill.
You should write a report that, as a minimum:
• discusses in detail the above tasks and
• contains well presented R code developed in responding to the tasks.
You are advised, but not obliged, to use R Markdown. We will certainly help you to get
started with R Markdown, if you ask. You can use LATEX from within R Markdown, but you
are under no obligation to use LATEX.
You should use clearly readable fonts and maintain sensible margins and line spacing through-
out. Your work should not look squashed. The default settings of R Markdown are completely
acceptable.
4.3 Coursework 2: What You Need to Submit
Your group needs to submit the following file electronically using the DLE:
• A Portable Document Format file MATH3609_Coursework_2_2020_2021.pdf containing
your report for Coursework 2.
All members of the group should take responsibility for seeing that your work
is submitted before the deadline.
It is also the responsibility of all members to check that the correct file has
been uploaded.
There have been occasions when R code has been submitted instead of the
required report. This error is easy to make, but will have severe consequences.
If anything is unclear, you should ask without delay.
39
学霸联盟
Coursework 1 and Coursework 2
Dr Nathan Broomhead and Dr Julian Stander
Module Taught in Semester 2, Academic Year 2020/21
1 General Information about Coursework 1 and
Coursework 2
Please read the following points before attempting the courseworks:
• The deadline for both Coursework 1 and Coursework 2 is 10 am on Thursday, 22nd
April, 2022. You should submit Coursework 1 and Coursework 2 separately through
the MATH3609 Optimization, Networks and Graphs DLE site. Your submission will be
marked anonymously.
• Coursework 1 and Coursework 2 are Group Courseworks. Please work in
a self-assigned group of up to four people. You will need to register your
group on the DLE in the usual way.
Within each coursework, each member of the group will receive the same mark, unless
any member chooses to make use of the Peer Assessment option. You should keep notes
of all your group meetings to use as evidence in case you choose to make use of the Peer
Assessment option. If you wish to make use of the Peer Assessment option, you will need
to make an appointment to see Dr Nathan Broomhead or Dr Julian Stander.
• You should give due consideration to your personal time management to ensure that
your work is submitted in plenty of time prior to the deadline. For further information,
please consult the Student Handbook
https://www.plymouth.ac.uk/your-university/governance/student-handbook,
https://www.plymouth.ac.uk/student-life/your-studies/essential-information/exams/exam-rules-and-
regulations/extenuating-circumstances,
https://www.plymouth.ac.uk/uploads/production/document/path/18/18487/Extenuating_Circumstances_
Policy_and_Procedures_20-21_v.3.pdf and
https://liveplymouthac.sharepoint.com/sites/x70/SitePages/Student-QA.aspx#q-a-session-1-safety-net-
video%E2%80%8B%E2%80%8B%E2%80%8B%E2%80%8B%E2%80%8B%E2%80%8B%E2%80%8B.
• Each member of the group can expect to spend up to around 15 hours on Coursework 1
and up to around 15 hours on Coursework 2. These suggested times do not include the
preparation work that you will need to do.
• Each coursework counts for 15% = 12 × 30% of your final mark on this module. Marks
will be assigned according to the marking grids on pages 26 and 30.
1
• Marked scripts will be returned within 20 working days of the submission date, and,
in particular, you will get full feedback on your work by Friday, 21st May, 2021, unless
extensions are used.
• You are reminded of the University’s Academic Regulations:
Academic offences occur when activity is undertaken which could confer an unfair
advantage to any candidate(s) in assessment. The University recognises the following
(including any attempt to carry out the actions described) as academic offences, regardless
of intent:
a. Plagiarism, which is copying or paraphrasing of other people’s work or ideas into a
submitted assessment without full acknowledgement. More information on plagiarism
is available here:
https://www.plymouth.ac.uk/student-life/your-studies/essential-information/regulations/plagiarism
b. Collusion, which is unauthorised collaboration of students (or others) in producing a
submitted assessment. The offence of collusion occurs if a student copies any part of
another student’s work, or allows their own work to be copied. Collusion also occurs
if other people contribute significantly to work that a student submits as their own.
The complete list of regulations can be found here:
https://www.plymouth.ac.uk/student-life/your-studies/essential-information/regulations
All group members automatically confirm that they have understood the University’s
policy on plagiarism and collusion when they submit each coursework.
We now state the relevant MATH3609 Optimization, Networks and Graphs Assessed Learning
Outcomes (ALOs) for this assignment.
At the end of the module the learner will be expected to be able to:
ALO1 Develop and analyse methods for the solution of optimisation problems;
ALO3 Demonstrate an understanding of the underlying concepts of graph theory;
ALO4 Use computer software for the solution of optimisation problems and report the
results.
You should keep these ALOs in mind when working on these courseworks.
2
2 Image Restoration: Preparation Work
You should study this section carefully. There has also been a lecture on it.
2.1 The Basic Idea
Let us begin by reading in an image and displaying it.
# Set the working directory
# This is where the image R_image_binary.csv is stored
#
# Please note that your working directory will be different
# from ours (only one of which is needed)
#
## setwd("/Users/julianstander/Dropbox/MATH3609_CW_2020_2021")
#
setwd("/Users/nbroomhead/Dropbox/MATH3609_CW_2020_2021")
#
R_logo <- as.matrix(read.table("R_image_binary.csv",
sep = ",")) # The image is thought of as a matrix
#
image(R_logo,
col = c("black", "white"))
0.0 0.2 0.4 0.6 0.8 1.0
0.
0
0.
2
0.
4
0.
6
0.
8
1.
0
Please note that this image is a binary image. It comprises two ‘colours’, black and white. In
Coursework 1 and Coursework 2 we work exclusively with binary images for simplicity.
We will discuss the axes labels later.
3
Unfortunately, often images are received in a degraded or noisy form such as
0.0 0.2 0.4 0.6 0.8 1.0
0.
0
0.
2
0.
4
0.
6
0.
8
1.
0
This could be because they are observed and then transmitted by a satellite and so are
subject to atmospheric distortions. Similarly, medical images are often distorted in some way;
just think of an ultrasound image of a baby or even an x-ray or magnetic resonance image
when there has been some patient movement.
4
The aim of image restoration is to try to restore the first image
0.0 0.2 0.4 0.6 0.8 1.0
0.
0
0.
2
0.
4
0.
6
0.
8
1.
0
from the degraded image
0.0 0.2 0.4 0.6 0.8 1.0
0.
0
0.
2
0.
4
0.
6
0.
8
1.
0
5
2.2 Some Mathematical Notation
Images are represented as matrices. Here is a small example to fix ideas:
x <- matrix(c(0, 0, 0,
0, 0, 1,
1, 1, 1,
0, 0, 1),
nrow = 4,
byrow = TRUE)
#
x
## [,1] [,2] [,3]
## [1,] 0 0 0
## [2,] 0 0 1
## [3,] 1 1 1
## [4,] 0 0 1
#
# Dimension
#
dim(x)
## [1] 4 3
nrow(x)
## [1] 4
ncol(x)
## [1] 3
#
# Specific
#
x[1, 2]
## [1] 0
x[3, 1]
## [1] 1
So the image that is stored in the R object x comprises 4 rows and 3 columns. The value
taken by the image in row 1, column 2 is 0 corresponding to black, while the value taken
by the image in row 3, column 1 is 1 corresponding to white. We say that the above image
comprises 4× 3 = 12 pixels. We shall use i to index the pixels: here i = 1, . . . , n = 12, where n
is the total number of pixels. The value taken at pixel i is denoted xi, where xi can take
the values 0 for black or 1 for white: xi ∈ {0,1}. This enables us to write the whole image as
the vector x = (x1, . . . , xn).
6
We label the pixels column by column:
Columns
R
ow
s
i = 1 or ( 1 , 1 )
i = 2 or ( 2 , 1 )
i = 3 or ( 3 , 1 )
i = 4 or ( 4 , 1 )
i = 5 or ( 1 , 2 )
i = 6 or ( 2 , 2 )
i = 7 or ( 3 , 2 )
i = 8 or ( 4 , 2 )
i = 9 or ( 1 , 3 )
i = 10 or ( 2 , 3 )
i = 11 or ( 3 , 3 )
i = 12 or ( 4 , 3 )
1 2 3
1
2
3
4
So that
x
## [,1] [,2] [,3]
## [1,] 0 0 0
## [2,] 0 0 1
## [3,] 1 1 1
## [4,] 0 0 1
is the vector
as.vector(x)
## [1] 0 0 1 0 0 0 1 0 0 1 1 1
or x = (x1, x2, . . . , x12) = (0,0,1,0,0,0,1,0,0,1,1,1).
7
Please note that in
Columns
R
ow
s
i = 1 or ( 1 , 1 )
i = 2 or ( 2 , 1 )
i = 3 or ( 3 , 1 )
i = 4 or ( 4 , 1 )
i = 5 or ( 1 , 2 )
i = 6 or ( 2 , 2 )
i = 7 or ( 3 , 2 )
i = 8 or ( 4 , 2 )
i = 9 or ( 1 , 3 )
i = 10 or ( 2 , 3 )
i = 11 or ( 3 , 3 )
i = 12 or ( 4 , 3 )
1 2 3
1
2
3
4
we also show the row and column numbers in brackets. We will not usually use these, but
they are added here for information. A point to notice is that the column labels increase
from left to right, while the row labels increase from top to bottom which is different from
standard Cartesian co-ordinates.
8
Let’s consider how R draws an image using its image function
x
## [,1] [,2] [,3]
## [1,] 0 0 0
## [2,] 0 0 1
## [3,] 1 1 1
## [4,] 0 0 1
#
image(x,
col = c("black", "white"))
0.0 0.2 0.4 0.6 0.8 1.0
−
0.
2
0.
2
0.
6
1.
0
#
# Or with better labels
#
n_row <- nrow(x) # Number of rows in x
n_row
## [1] 4
#
n_col <- ncol(x) # Number of columns in x
n_col
## [1] 3
#
image(1:n_row, # Label the rows
1:n_col, # Label the columns
x,
col = c("black", "white"),
xlab = "Row number",
9
ylab = "Column number",
axes = FALSE)
#
axis(1, 1:n_row)
axis(2, 1:n_col)
Row number
Co
lu
m
n
nu
m
be
r
1 2 3 4
1
2
3
So R thinks of the horizontal axis as representing row number, and the vertical axis as
representing column number, just like standard Cartesian co-ordinates: (row i, column j) =(horizontal co-ordinate x,vertical co-ordinate y). This means that the first column of the
matrix x is plotted horizontally at the bottom, the middle column of x is plotted horizontally
in the middle, and the final column of x is plotted horizontally at the top.
These plotting details are not important for what will come. However, we give them here for
information.
10
2.3 A Little More Mathematical Notation
We have seen that in practice we observe a degraded version of a ‘true’ image such as:
0.0 0.2 0.4 0.6 0.8 1.0
0.
0
0.
2
0.
4
0.
6
0.
8
1.
0
We refer to this observed or degraded image using the vector y = (y1, . . . , yn), where n is the
number of pixels. Here yi ∈ {0,1}, i = 1, . . . , n.
11
2.4 Some Literature and Other Irrelevant Information
In scientific writing, we should always give credit where credit is due. Many people have
contributed to statistical image restoration. The following are but two references among a
great many:
Besag, J. E. (1986) On the statistical analysis of dirty pictures (with discussion). Journal of
the Royal Statistical Society B, 48, 259–302.
Geman, D. and Geman, S. (1984) Stochastic relaxation, Gibbs distributions and the Bayesian
restoration of images. IEEE Transactions on Pattern Analysis and Machine Intelligence, 6,
721–41.
Coursework 1 and Coursework 2 are based on the work of
Greig, D. M., Porteous, B. T. and Seheult, A. H. (1989) Maximum a posteriori estimation
for binary images. Journal of the Royal Statistical Society B, 51, 271–279.
Greig, Porteous and Seheult (1989) assume black is coded as 1 and white as 0. In Coursework
1 and Coursework 2, we use take 1 to represent white and 0 to represent black, as this is
more natural and usual.
The work of Greig, Porteous and Seheult (1989) was expanded on in these works:
Jubb, M. D. (1989) Image Reconstruction. PhD thesis. University of Bath.
Stander, J. (1992) Some Topics in Statistical Image Analysis. PhD thesis. University of Bath.
Stander, J. and Silverman, B. W. (1994) Temperature schedules for simulated annealing.
Statistics and Computing, 4, 21–32.
Franconi, L. and Jennison, C. (1997) Comparison of a genetric algorithm and simulated
annealing in an application to statistical image reconstruction. Statistics and Computing, 7,
193–207.
No knowledge of any of this literature is assumed in either Coursework 1 or
Coursework 2.
We end this section by remarking that images are not usually binary. Often each pixel can
take a grey level so that xi ∈ {0,1, . . . ,255} where 0 is black and 255 is white for example,
or a colour, which is more complicated to represent. We work with binary images to keep
life relatively simple, as this may be the first time that you have met image reconstruction
problems.
12
2.5 Statistical Models
Image restoration proceeds by specifying statistical models for possible true images x =(x1, . . . , xn) and for degraded images y = (y1, . . . , yn).
Please note that this is not a statistics module, so not all the details that we present here are
needed for you to be able to do Coursework 1 or Coursework 2. We include these details
to explain how the image restoration problem can be expressed as an optimization
problem. It is with that optimization problem that we shall work.
2.6 A Statistical Model for y given x
We will assume the following statistical model for yi given xi:
Pr(yi = 0 ∣xi = 0) = 1 − p (no error made, no swap)
Pr(yi = 1 ∣xi = 0) = p (error made, swap from 0 to 1)
Pr(yi = 0 ∣xi = 1) = p (error made, swap from 1 to 0)
Pr(yi = 1 ∣xi = 1) = 1 − p (no error made, no swap)
Essentially, we are assuming that the colour at pixel i is swapped from black to white, or
from white to black, with probability p ∈ [0,1]; otherwise, it is not changed. For simplicity,
we assume that p is known.
We assume that the error mechanism at each pixel is independent of the error mechanisms at
every other pixel. This enables us to write down Pr(y ∣x):
Pr(y ∣x) = n∏
i=1 Pr(yi ∣xi).
Let’s recall that the degraded image y is observed, while the ‘true’ image x is unknown and
needs to be estimated. Thought of as a function of the unknown image vector x, Pr(y ∣x)
is the likelihood function that you may have met in MATH2602 Statistical Inference and
Regression or MATH3613 Data Modelling. We shall write L(y ∣x) = Pr(y ∣x).
13
Now let’s write f(yi ∣xi) = Pr(yi ∣xi), so that f(1 ∣0) = Pr(yi = 1 ∣xi = 0) for example. We
have
L(y ∣x) = Pr(y ∣x)
= n∏
i=1 Pr(yi ∣xi)
= n∏
i=1 f(yi ∣1)xif(yi ∣0)1−xi ,
since f(yi ∣1)xif(yi ∣0)1−xi = f(yi ∣1) if xi = 1 and f(yi ∣0) if xi = 0.
From this it follows that
L(y ∣x) = n∏
i=1 f(yi ∣0) × {f(yi ∣1)f(yi ∣0)}
xi
. (1)
Later, we will use λi, which is defined as follows:
λi = log{f(yi ∣1)
f(yi ∣0)} . (2)
Please note that λi only depends on the observed yi and not on xi. Please also note that, if
yi = 0, then
λi = log{f(yi ∣1)
f(yi ∣0)} = log{f(0 ∣1)f(0 ∣0)} = log ( p1 − p) ,
while, if yi = 1, then
λi = log{f(yi ∣1)
f(yi ∣0)} = log{f(1 ∣1)f(1 ∣0)} = log (1 − pp ) .
14
2.7 A Statistical Model for x
To define a statistical model for x, we need to explain what it means for two pixels to be
neighbours. This can be done by the following picture
Columns
R
ow
s
i = 1 or ( 1 , 1 )
i = 2 or ( 2 , 1 )
i = 3 or ( 3 , 1 )
i = 4 or ( 4 , 1 )
i = 5 or ( 1 , 2 )
i = 6 or ( 2 , 2 )
i = 7 or ( 3 , 2 )
i = 8 or ( 4 , 2 )
i = 9 or ( 1 , 3 )
i = 10 or ( 2 , 3 )
i = 11 or ( 3 , 3 )
i = 12 or ( 4 , 3 )
1 2 3
1
2
3
4
The red and blue lines indicate neighbouring pixles; we use blue to indicate horizontal
neighbours and red to indicate vertical neighbours, although this distinction is not important.
As an example we see that the neighbours of pixel i = 6 are 5, 10, 7 and 2. It turns out
that it simpler to assume that the image is wrapped on a torus, so that, for example the
neighbours of pixel i = 1 are 4 (wrapped top to bottom), 5, 2 and 9 (wrapped left to right).
This assumption has very little practical effect for large images, but it does mean that every
pixel has four neighbours.
15
These neighbourhood relationships lead to a symmetric adjacency matrix A, which here
will be of dimension n × n = 12 × 12:
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## [1,] 0 1 0 1 1 0 0 0 1 0 0 0
## [2,] 1 0 1 0 0 1 0 0 0 1 0 0
## [3,] 0 1 0 1 0 0 1 0 0 0 1 0
## [4,] 1 0 1 0 0 0 0 1 0 0 0 1
## [5,] 1 0 0 0 0 1 0 1 1 0 0 0
## [6,] 0 1 0 0 1 0 1 0 0 1 0 0
## [7,] 0 0 1 0 0 1 0 1 0 0 1 0
## [8,] 0 0 0 1 1 0 1 0 0 0 0 1
## [9,] 1 0 0 0 1 0 0 0 0 1 0 1
## [10,] 0 1 0 0 0 1 0 0 1 0 1 0
## [11,] 0 0 1 0 0 0 1 0 0 1 0 1
## [12,] 0 0 0 1 0 0 0 1 1 0 1 0
Do you understand why?
16
We can now define a probability mass function on images x as
Pr(x) = exp (−β × number of discrepant neighbouring pairs in x)∑images x′ exp (−β × number of discrepant neighbouring pairs in x′) ,
in which β ≥ 0 is a constant that we will discuss later. The sum over all images x′ comprises
2n terms, where n is the number of pixels, and can be very hard to compute when n is large.
As this sum does not depend on the argument x of Pr(x), we can discard it, if we replace
the equality with a ‘proportional to’ sign. So we write
Pr(x)∝ exp (−β × number of discrepant neighbouring pairs in x) .
For the image x that we have been working with:
## [,1] [,2] [,3]
## [1,] 0 0 0
## [2,] 0 0 1
## [3,] 1 1 1
## [4,] 0 0 1
the number of discrepant neighbouring pairs is 10:
Columns
R
ow
s
0
0
1
0
0
0
1
0
0
1
1
1
1 2 3
1
2
3
4
Hence, the probability of this image is Pr(x) = C exp(−β × 10), where C is a constant.
17
Similarly, the number of discrepant neighbouring pixels in the all black and the all white
images is, unsurprisingly, 0:
Columns
R
ow
s
0
0
0
0
0
0
0
0
0
0
0
0
1 2 3
1
2
3
4
Columns
R
ow
s
1
1
1
1
1
1
1
1
1
1
1
1
1 2 3
1
2
3
4
This means that the all black and the all white images are the most probable with probabilities
C exp(β × 0) = C each.
18
The chess-board image is the least probable, as it has the highest number of discrepant
neighbouring pixels:
Columns
R
ow
s
1
0
1
0
0
1
0
1
1
0
1
0
1 2 3
1
2
3
4
There are 20 discrepant pairs, meaning that the probability of the chess board image is
C exp(−β × 20).
You may be asking: what does β ≥ 0 do? The larger the value of β, the higher the relative
probabity given to the all black and to the all white image. In other words, higher β values
encourage ‘smoother’ images. Because of this, sometimes β is referred to as a smoothing
parameter.
Finally, we remark that our probability model for images x,
Pr(x)∝ exp (−β × number of discrepant neighbouring pairs in x) ,
is essentially the same as the Ising model used a lot in physics and statistical mechanics
to represent substances such as magnets. To quote Federico Camia: ‘The Ising model,
introduced by Lenz in 1920 to describe ferromagnetism, is one of the most studied models
of statistical mechanics. Its two dimensional version has played a special role in rigorous
statistical mechanics since Peierls’ proof of a phase transition in 1936 and Onsager’s derivation
of the free energy in 1944.’
19
It can be shown that the
number of discrepant neighbouring pairs in x = 12 n∑i=1 n∑j=1aij (xi − xj)2 ,
in which x = (x1, . . . , xn) is an image and aij is the (i, j)-th element of the adjacence matrix
A such that
aij = {1 if pixels i and j are neighbours0 if pixels i and j are not neighbours.
A consequence of this result is that
Pr(x) ∝ exp (−β × number of discrepant neighbouring pairs in x)
= exp(−β × 12 n∑i=1 n∑j=1aij (xi − xj)2)
= exp(−β2 n∑i=1 n∑j=1aij (xi − xj)2) . (3)
Another irrelevant fact is that more general probability mass functions can be defined by
replacing (xi − xj)2 with φ (xi − xj) where the function φ satisfies certain properties.
20
2.8 A Statistical Model for x given y
It makes sense to consider Pr(x ∣ y), as we are given the observed image y. Pr(x ∣ y) can be
found using Bayes’ theorem as follows:
Pr(x ∣ y) ∝ L(y ∣x) ×Pr(x)
∝ n∏
i=1 f(yi ∣0) × {f(yi ∣1)f(yi ∣0)}
xi ×Pr(x) from (1)
∝ n∏
i=1 f(yi ∣0) × {f(yi ∣1)f(yi ∣0)}
xi × exp(−β2 n∑i=1 n∑j=1aij (xi − xj)2) from (3) (4)
In the language of Bayesian statistical inference, Pr(x) is the prior probability mass
function of all images x before observing the degraded image y. Pr(x) expresses our belief
about x before seeing the data in terms of probabilities. Pr(x ∣ y) is our belief about x after
seeing the data y and is referred to as our posterior probability mass function.
Please note that, as mentioned, it does not matter if you do not follow all these details,
although they are not hard. The important part is that image restoration leads to an
optimization problem, as we shall see.
21
2.9 The Restoration of a Binary Image as an Optimization Problem
The restored image is defined as the most probable image under the posterior Pr(x ∣ y).
That is, it is the image x that maximizes Pr(x ∣ y). For this reason it is referred to as the
maximum a posteriori or MAP image.
From (4) we see that the restored image maximizes
Pr(x ∣ y)∝ n∏
i=1 f(yi ∣0) × {f(yi ∣1)f(yi ∣0)}
xi × exp(−β2 n∑i=1 n∑j=1aij (xi − xj)2)
over all images x.
We know that it’s very often easier to work with the logarithm of a function, rather than
the function itself, when maximizing. This is because the logarithm of a function and the
function itself are maximized at the same point, and also because taking logarithms turns
products into sums, which are easier to deal with. So, the restored image maximizes
log {Pr(x ∣ y)} = additive constant D + n∑
i=1 log [f(yi ∣0) × {f(yi ∣1)f(yi ∣0)}
xi] − β2 n∑i=1 n∑j=1aij (xi − xj)2
= D + n∑
i=1 log {f(yi ∣0)} + n∑i=1 log [{f(yi ∣1)f(yi ∣0)}
xi] − β2 n∑i=1 n∑j=1aij (xi − xj)2
= D + n∑
i=1 log {f(yi ∣0)} + n∑i=1 xi log{f(yi ∣1)f(yi ∣0)} − β2 n∑i=1 n∑j=1aij (xi − xj)2
= D + n∑
i=1 log {f(yi ∣0)} + n∑i=1 xiλi − β2 n∑i=1 n∑j=1aij (xi − xj)2 (5)
by the definition of λi in (2).
Please note that ∑ni=1 log {f(yi ∣0)} does not depend on x, and so can be considered a constant.
It therefore follows that we need to maximize
n∑
i=1 xiλi − β2 n∑i=1 n∑j=1aij (xi − xj)2
over x = (x1, . . . , xn), or, equivalently, to minimize
− n∑
i=1 xiλi + β2 n∑i=1 n∑j=1aij (xi − xj)2
over all images x.
22
So, we have finally arrived at our image restoration optimization problem. Given a degraded
image y = (y1, . . . , yn), we can find our restored image by solving the problem:
minimize: − n∑
i=1 xiλi + β2 n∑i=1 n∑j=1aij (xi − xj)2 over images x = (x1, . . . , xn). (6)
It can be shown that the minimiation problem given by (6) is equivalent to the problem
minimize:
n∑
i=1(1 − xi)max(0, λi) + n∑i=1 ximax(0,−λi) + β2 n∑i=1 n∑j=1aij (xi − xj)2 (7)
over images x, in which max(a, b) means the maximum of a and b, so that, for example,
max(0,1) = 1 and max(0,−1) = 0.
This quantity is never negative.
23
2.10 Solving the Minimization Problem using the Ford-Fulkerson
Algorithm
We will now see how Greig, Porteous and Seheult (1989) solved the minimization problem
(7), namely,
minimize:
n∑
i=1(1 − xi)max(0, λi) + n∑i=1 ximax(0,−λi) + β2 n∑i=1 n∑j=1aij (xi − xj)2
over images x, using the Ford-Fulkerson algorithm. To do this, Greig, Porteous and Seheult
(1989) defined a network with capacities as follows:
2.10.1 Definition of the Network with Capacities
• There are n + 2 vertices:
– A source s;
– A sink t;
– A vertex for each of the n pixels, indexed by i or j below.
• If λi > 0, there is an arc (directed edge) (s, i) from the source s to i with capactiy
csi = λi; otherwise
• If λi < 0, there is an arc (i, t) from i to the sink t with capactiy cit = −λi;
• There is an undirected edge {i, j} between two internal vertices i and j (that corre-
spond to pixels i and j) with capacities cij = β aij. This can be thought of as two arcs,
one between i and j, and the other between j and i, both of which have capacities β aij .
• There are no other arcs.
24
2.10.2 Definition of a Cut
Let us say that x = (x1, . . . , xn) is an image that is being proposed as a reconstruction of the
observed y = (y1, . . . , yn).
• Partitioning the Network
We can use x to partition the network into two sets:
W (x) = {s} ∪ {i ∶ xi = 1}, source s and all white pixels
B(x) = {i ∶ xi = 0} ∪ {t}, all black pixels and the sink t.
This partition of the vertices in the network into the sets W (x) and B(x) defines a cut.
• The Cut Capacity
The associated cut capacity is
C(x) ∶= C(∂+(W (x))) = ∑
kl∈∂+(W (x)) ckl. (8)
Notice that the arcs of ∂+(W (x)) are of three types: arcs starting at s; arcs ending at t; arcs
between intermediate vertices.
The cut capacity (8) can also be written as
C(x) = ∑
k∈W (x) ∑l∈B(x) ckl.
It turns out that
C(x) = n∑
i=1(1 − xi)max(0, λi) + n∑i=1 ximax(0,−λi) + β2 n∑i=1 n∑j=1aij (xi − xj)2 ,
which is exactly the function that we wanted to minimize in (7) to restore the image from y.,
So, the image restoration problem becomes a max-flow/min-cut problem: we find the cut
that minimizes the cut capacity C(x) and this determines x. We can do this using the
Ford-Fulkerson algorithm (or modern improvements of it). ,
We will expand on and explore the above in Coursework 1 and Coursework 2.
25
3 Coursework 1
Marks will be assigned for Coursework 1 according to this marking grid.
MATH3609 Optimization, Networks and Graphs:
Coursework 1 Marking Grid
Assessment Area Maximum Awarded Feedback
Mark Mark
Task 1: Correctness and va-
lidity of argument. Read-
ability of explanation.
15
Task 2: Correctness and va-
lidity of argument. Read-
ability of explanation.
20
Task 3: Correctness and va-
lidity of argument. Read-
ability of explanation.
10
Task4: Correctness and va-
lidity of argument. Read-
ability of explanation.
25
Task 5: Correctness and va-
lidity of argument. Read-
ability of explanation.
15
Task 6: Correctness and va-
lidity of argument. Read-
ability of explanation.
15
Total1 100
1This mark is provisional. Like all marks it is considered by an External Examiner and an Assessment Panel.
26
3.1 Coursework 1 Tasks
Coursework 1 comprises a number of tasks. You need to produce a document containing
your solutions to these tasks. These may be hand-written and then scanned or typed up
using, for example, LATEX. It is your responsibility to make sure that they are readable.
The page limit is fifteen pages. Please note that this is the limit and should not be
considered as the target. Well written, clear and concise mathematics will receive high credit.
Please arrange your report under subheadings referring to the tasks, that is, Task 1,
Task 2, . . .. This will help you to respond systematically to all parts of Coursework 1.
The marks for each task are shown in the grid on page 26.
Task 1: Let G = (V,E) be a graph with vertex set V = {1,2, . . . , n} and let A be the
adjacency matrix of G with (i, j)-th entry aij. Suppose x∶V → {0,1} is a function which
assigns to each vertex i a number xi, which is either 0 or 1. A pair of neighbouring vertices i
and j are called discrepant if xi ≠ xj.
Explain why the number of discrepant neighbouring (unordered) pairs is given by the formula:
1
2
n∑
i=1
n∑
j=1aij (xi − xj)2 .
Task 2: Explain why the minimization problem given by (6) on page 23 is equivalent to the
problem
minimize:
n∑
i=1(1 − xi)max(0, λi) + n∑i=1 ximax(0,−λi) + β2 n∑i=1 n∑j=1aij (xi − xj)2
over images x, in which max(a, b) means the maximum of a and b, so that, for example,
max(0,1) = 1 and max(0,−1) = 0.
Task 3: Prove that
n∑
i=1(1 − xi)max(0, λi) + n∑i=1 ximax(0,−λi) + β2 n∑i=1 n∑j=1aij (xi − xj)2
is never negative.
Task 4: Prove that
C(∂+(W (x))) = n∑
i=1(1 − xi)max(0, λi) + n∑i=1 ximax(0,−λi) + β2 n∑i=1 n∑j=1aij (xi − xj)2 .
Hint: The arcs of ∂+(W (x)) can be split into three types, namely, arcs that start at s, arcs
that end at t and arcs between intermediate vertices. Consider the sums of the capacities
arcs of each type separately.
27
Task 5: Show that every cut K of the network can be written as
K = ∂+(W (x))
for some image x, and deduce using this and the previous tasks that the Ford-Fulkerson
algorithm can be used to find a restored image.
Task 6: Is there always a unique solution x to the minimization problem (i.e. a unique
restored image)? Illustrate your answer with an example.
28
3.2 Coursework 1: What You Need to Submit
Your group needs to submit the following file electronically using the DLE:
• A Portable Document Format file MATH3609_Coursework_1_2020_2021.pdf containing
your report for Coursework 1.
All members of the group should take responsibility for seeing that your work
is submitted before the deadline.
It is also the responsibility of all members to check that the correct file has
been uploaded and that it is readable.
If anything is unclear, you should ask without delay.
29
4 Coursework 2
Marks will be assigned for Coursework 2 according to this marking grid.
MATH3609 Optimization, Networks and Graphs:
Coursework 2 Marking Grid
Assessment Area Maximum Awarded Feedback
Mark Mark
Image Restoration
Tasks: Defining the set-
up. Correctness, clarity
and presentation of code.
25
Image Restoration
Tasks: Optimizing,
presenting the results,
choosing β. Depth of
understanding and quality
of description of insights
gained.
50
Report: display results
clearly, concisely, system-
atically and professionally,
using correct spelling and
grammar. This is important
for employability.
25
Total2 100
2This mark is provisional. Like all marks it is considered by an External Examiner and an Assessment Panel.
30
4.1 Coursework 2 Tasks
Coursework 2 comprises a number of tasks. You need to produce a report of your work
following the instructions below.
Your report should contain well presented and annotated R code.
The page limit is twenty pages, including R code and figures. Please see page 39
for more details. Please do not submit an additional appendix as it will not be considered.
Reports that contain irrelevant or uninteresting discussion or R code will be penalized. Please
note that twenty pages is the limit and should not be considered as the target. Well written,
concise reports will receive high credit.
You may also omit lengthy, technical R code, if you feel that it does not add to
the argument that you are making. Essential code should, however, be shown.
Please arrange your report under subheadings referring to the tasks, that is, Task 1,
Task 2, . . .. This will help you to respond systematically to all parts of Coursework 2.
31
4.1.1 Implementing the Ford-Fulkerson Algorithm to Restore the
Image in R: Getting Started
Following on from the information in the section ‘Image Restoration: Preparation Work’ we
will now produce an image restoration from observed data y using the Ford-Fulkerson
algorithm implemented in R.
First, we need to read y into R. The data that we will use are supplied in the file
R_degraded_image_0.2.csv. We can read in and display these data in the same was
as we did at the beginning of this document. Here is the code again to make sure that you
get started properly:
y <- as.matrix(read.table("R_degraded_image_0.2.csv",
sep = ","))
#
image(y,
col = c("black", "white"))
0.0 0.2 0.4 0.6 0.8 1.0
0.
0
0.
2
0.
4
0.
6
0.
8
1.
0
#
n_row <- nrow(y)
n_col <- ncol(y)
#
c(n_row, n_col)
## [1] 23 28
#
y_vec <- as.vector(y)
#
length(y_vec)
## [1] 644
32
We can label the axes in a possibly nicer way:
image(1:n_row,
1:n_col,
y,
xlab = "Row number",
ylab = "Column number",
main = "Data: degraded image y",
col = c("black", "white"))
5 10 15 20
5
10
15
20
25
Data: degraded image y
Row number
Co
lu
m
n
nu
m
be
r
• It is not hard to imagine that y is a degraded version of the image in R_logo seen above.
• It is known that p = 0.2. (In practice, p and any other unknown parameter or parameters
can also be estimated form the data by assigning prior distributions to them and sampling
from the resulting posterior using Markov chain Monte Carlo methods as discussed in
MATH3613 Data Modelling. This is beyond the scope of Coursework 2.)
• We will assume for now that β = 1.
33
4.1.2 Defining the Vector of λ Values
Task 1:
Using the function ifelse or otherwise, define the vector of λ values (λ1, . . . , λn) in R and
show a few λ values. Identify and record in R which λi values are positive and which λi values
are negative using the which function, storing the results in vectors which_lambda_positive
and which_lambda_negative respectively.
4.1.3 Defining the Network with Capacities
We need to define a data frame comprising three variables or columns called from, to and
capacity. This data frame defines the capacities in the network. If, for example, we wish to
specify ckl, then the variable from would take the value k, the variable to would take the
value l and the variable capacity would take the value ckl.
One way to define the required data frame of capacities is to construct three data
frames df_source, df_sink and df_internal, and to combine df_source, df_sink and
df_internal by binding their rows into one data frame df using rbind.
The data frame df_source records the capacities of arcs coming from the source s.
Task 2: Define df_source in R. The square brackets [] may help you here. In order to
check this, the first few rows of the result are given
## from to capacity
## 1 source 2 1.386294
## 2 source 3 1.386294
## 3 source 4 1.386294
## 4 source 6 1.386294
## 5 source 7 1.386294
## 6 source 8 1.386294
The data frame df_sink records the capacities of arcs going to the sink t.
Task 3: Define df_sink in R and show a part of it.
34
Before we define df_internal, we need to obtain the adjacency matrix A. This can be
done by using the function get_neighbourhood_matrix with which you are supplied in the
file of helper functions helper_functions.R; you do not need to understand the details of
these helper functions. Please note that get_neighbourhood_matrix is defined using other
functions, so that you’ll have to make sure that all the functions in the file of helper functions
helper_functions.R are loaded into R. We also need to extract the rows and columns of A
corresonding to neighbouring pixels:
A <- get_neighbourhood_matrix(y, n_row)
#
A[1:5, 1:5]
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0 1 0 0 0
## [2,] 1 0 1 0 0
## [3,] 0 1 0 1 0
## [4,] 0 0 1 0 1
## [5,] 0 0 0 1 0
#
neighbouring_pixels <- which(A == 1,
arr.ind = TRUE)
#
head(neighbouring_pixels)
## row col
## [1,] 2 1
## [2,] 23 1
## [3,] 24 1
## [4,] 622 1
## [5,] 1 2
## [6,] 3 2
35
The data frame df_internal records cij where i and j correspond to pixels i and j.
Task 4: Define df_internal in R and show a part of it. Define the data frame df using
rbind. Use the function graph_from_data_frame from the igraph package to create an
igraph graph called g from df. Produce a plot of your capacitated network, as shown here;
please note, that it’s very big (but rather nice)! Please note that your plot may be displayed in
a different way from the one here. You may find https://igraph.org/r/doc/plot.common.html
helpful.
source
2
23
24
622
1
3
25
623
4
26
624
5
27
625
6
28
626
7
29
627
8
30
628
9
31
629
10
32
630
11
33
631
12
34
632
13
35633
14 36
634
15
37
635
16
38
636
17
39 637
18
40
638
19
41
639
20
42
640
21
43
64122
44
642
45
643
46
64447
48
49
50
51
52
53
54
55
56
57
58
5960
6162
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
8283
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124 125
126
127
128
129
130
131132
133
134
135
136
137
138
139
140
141
142
143
144
145
146 147
148
149
150
151
152153
154
155
156
157
158
159
160
161 162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195196
197198
199
200
201
202
203
204
205
206
207 208
209
210
211
212
213
214
215
216
217218219220
221
222
223
224
225
226
227
228
229
230 231
232 233
234
235
236
237
238
239
240
241242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264265
266
267
268
269
270
271
272
273
274 275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300 301
302
303
304
305
306
307
308
309
310
311
312313314
315
316
317
318
319
320
321 322
323
324
325
326
327
328
329330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384385
386
387
388
389
390
391
392
393
394
395
396
397
398
399400
401
402
40
404
405
406
407
408409
410
411
412
413
414
415
416
417
418
419
420
421
422423
424
425
426
427
428
429
430
431
432
433
434
435
436 437
438
439
440
441
442
443
444
445446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
47
472
473
474
475
476
477
478
79
480
481
482
483
484
485486487
488
489
490
491 492
493
494
495
496
497498
499
500
501
502
503
504
505
506
507
508
509
510
511
512513
514
515
516
517
518
519
520521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542543
544
545
546
547
548
549550
551
552553
554
555
556
557
558
559
560
561
562
563
564
565566567
568
569
570
571
572573574
575
576
577
578
579
580
581
582
583
584
585
586 587
588589590
591
592
593
594
595
596
597
598
599
600601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
sink
36
4.1.4 Producing the Image Restoration by Maximizing Flow or Minimizing Cut Capacity
The function max_flow provided by the igraph package finds the maximum flow in a graph
or equivalently the minimum cut capacity. max_flow effectively implements an enhanced
modern version of the Ford-Fulkerson algorithm. Amongst other output, max_flow
returns a partition of the vertices into two sets that it calls partition1 and partition2.
partition1 is the set W (x) and partition2 is the set B(x) defined above, corresponding to
the ‘maximum flow’/‘minimum cut’ x. partition1 and partition2 are quite complicated
R objects. However, the part of interest to us can be extracted using the names function;
a character vector is given by names. Once names has been applied and any unnecessary
element has been removed, the resulting character vector can be turned into pixels numbers
using the as.numeric function.
Task 5: Use max_flow to maximize the flow in g. Report the value of the maximum flow.
Working with partition1 find the numerical values of the white pixels in the set W (x),
corresponding to the max-flow/min-cut, and report the first six of these. Similarly, find the
numerical values of the black pixels in the set B(x) and report the first six.
Task 6: Use these pixel values to define the reconstructed image. An easy way (but not the
only way) to do this is to use rep to define a vector, x_map_vec say, of the same length as
y_vec comprising NAs, and then to assign the value 0 to pixels that have been identified as
black and the value 1 to pixels that have been identified as white. Then turn the vector into
a matrix, remembering that pixels are labelled column by column; please see the help file for
matrix, if necessary. Display your image. Here is the one that we obtained.
0.0 0.2 0.4 0.6 0.8 1.0
0.
0
0.
2
0.
4
0.
6
0.
8
1.
0
37
4.1.5 Choosing β
This reconstruction x obtained from the observed image y was obtained using β = 1. In fact,
the choice of β is in the hands of the user.
Task 7:
Discuss and illustrate how you could choose β.
There is no right answer here. What we are looking for is sensible suggestions for
the choice of β, discussed clearly and critically. Your discussion should be illus-
trated with useful figures, as appropriate. Concise answers comprising sensible
and well illustrated suggestions will receive high credit.
38
4.2 Coursework 2: Report Production
The ability to report technical results clearly and concisely is an important skill.
You should write a report that, as a minimum:
• discusses in detail the above tasks and
• contains well presented R code developed in responding to the tasks.
You are advised, but not obliged, to use R Markdown. We will certainly help you to get
started with R Markdown, if you ask. You can use LATEX from within R Markdown, but you
are under no obligation to use LATEX.
You should use clearly readable fonts and maintain sensible margins and line spacing through-
out. Your work should not look squashed. The default settings of R Markdown are completely
acceptable.
4.3 Coursework 2: What You Need to Submit
Your group needs to submit the following file electronically using the DLE:
• A Portable Document Format file MATH3609_Coursework_2_2020_2021.pdf containing
your report for Coursework 2.
All members of the group should take responsibility for seeing that your work
is submitted before the deadline.
It is also the responsibility of all members to check that the correct file has
been uploaded.
There have been occasions when R code has been submitted instead of the
required report. This error is easy to make, but will have severe consequences.
If anything is unclear, you should ask without delay.
39
学霸联盟
