Hello, if you have any need, please feel free to consult us, this is my wechat: wx91due

COMP0210 2024/25: Research Computing with C++ Assignment 1

Errata

If you believe you have found an error in this assignment, please contact the module leader as soon as possible.

• Corrected signs in description of algorithm in section 0.1.3; you should subtract the changes to the image vector rather than add them.

• Correct test image named to TestSphere.pgm in section 4.4

• Corrected order of axes in FFTW library calls

Late Submission Policy

Please be aware of the UCL Late Submission Policy (Section 3.12).

Extenuating Circumstances

The course tutors are not allowed to grant Extenuating Circumstances for late submission. You must apply by submitting an Extenuating Circumstances form which can be found on the Extenuating Circumstances information website. The Department EC Panel will then decide whether or not to grant extenuating circumstances and email the module lead confirming that Extenuating Circumstances have been granted and providing a revised submission date. 

Collaboration and ChatGPT / AI Tools

This assignment must be you own work. You are not permitted to share or collaborate on code, nor are you permitted to use AI tools such as ChatGPT or github co-pilot to produce any code or question responses.

You may discuss the background of the assignment to help you understand the problem at hand, but not any solutions.

You are strongly encouraged to use the Moodle forums as a means of asking and responding to questions. You can also contact the lecturer directly via email or booking office hours.

Assignment Starter Download and Submission

The starter repository must be downloaded from the github link below, which will generate your personal github repository. You can use GitHub for your work if you wish, although we will only mark the code you submit on Moodle. Due to the need to avoid plagiarism, do not use a public GitHub repository for your work. We have provided a private repo unique to you accessible through Github Classroom. To access the starting project, using the following link:

https://classroom.github.com/a/LqDHSaBT.

After you accept the permissions, you will have access to a repository named

image-reconstruction-<gh_username>

for your particular github username. This will contain the starting project for your own work.

Important: It is vital you use the above link to access the starting project.

For this coursework assignment, you must submit by uploading a single Zip format archive to Moodle. You must use only the zip tool, not any other archiver, such as .tgz or .rar. Inside your zip archive, you must have a single top-level folder, whose folder name is your student number, so that on running unzip this folder appears. This top-level folder must contain all the parts of your solution. Inside the top-level folder should be a git repository named exactly COMP0210Assignment1. All code, images, results and documentation should be inside this git repository. You must include the .git folder otherwise we cannot mark your git commit log.

Your folder structure should look similar to (but not exactly the same as):

- 314159 // YOUR STUDENT NUMBER

- COMP0210Assignment1

- .git

- README.md

- source

- test

- ...

• Please remove the build folder before zipping and submitting.

• In order to have this folder structure I recommend starting by creating a folder with your student number, and then moving into that folder to clone the repository.

Important for cloning this repository: - You should clone this repository into a folder called COMP0210Assignment1 rather than the default name. This is so that everyone has the same folder name and the automated tools can access the folder and build your code.

You should therefore clone with the command:

git clone <your private github link> COMP0210Assignment1

N.B. We have provided a .devcontainer folder in this repository. To develop in the container, make sure that you have Docker running and then select “Reopen in Container” with the top level folder COMP0210Assignment1 open in VSCode.

You should test that your zip file produces this structure; if you are working on Windows, you should zip/unzip the files inside WSL, as the native Windows tool will produce a different folder structure when unzipped. Be sure to test your submission by unzipping the .zip file and checking the folder structure before you submit; will guarantee that your zip file will extract correctly during the submission process. You can use the command line tool zip in WSL / Linux terminals.

You will notice that there are multiple source folders, Optimisation (library code), app (appli cation code), and test (testing code). You should take care to choose an appropriate location for each piece of code you write. You can, and should, add new source and header files where appropriate. You should consider where they fit into the overall structure of the repository and abide by good practices regarding the grouping of related code in files, and the splitting of declarations (header files) and implementations (source files).

If writing tests involving paths (e.g. when testing loading a grid from a specific file) the test executables will be run from the top level directory. You may hard-code relative paths,provided they execute  correctly when run from the top level directory.

We recommend you use Visual Studio Code IDE with either Ubuntu, WSL running Ubuntu, or a devcontainer running Ubuntu 22.04 docker image for development (the correct Dockerfile is included for you in the repository). The devcontainer is recommended for Mac users to minimise compatibility issues. This is not a requirement but if you do use a different setup then you need to ensure that your code compiles and runs on Ubuntu 20.04 or 22.04, with g++ (enforcing C++17) and CMake 3.21 or above as this is the environment the markers will use to test the code. Your code should not rely on any external dependencies other than Catch2 for testing, and FFTW3 for the Fourier transforms.

This assignment requires installation of the Catch2 and FFTW3 libraries. I also strongly recommend having Valgrind or a similar tool available. All of these are already installed in the devcontainer, so consider using Docker with the container if you have any issues or want to avoid installs on your own machine.

Marks will be deducted for code that does not compile or run. Marks will also be deducted if code is poorly formatted. You can choose to follow your own formatting style but it must be applied consistently. See the Google C++ style guide for a good set of conventions regarding code formatting.

Naming and Build Structure

1. It is important to use the exact given names for any executable files and to use the correct command line interface described in the assignment, as these will be compiled and run automatically. You are free to add .cpp and .h files throughout the project, and name them yourself, but the executables must retain their names and interfaces.

2. All of your tests should be compiled into a single test executable, with the given name TestOptimisation, as in the initial CMake files. You may use as many additional .cpp files as you like, but they should still all compile down to a single executable with this name.

3. All executables must be placed in build/bin, as in the initial CMake files.

4. Where names of functions or classes etc. are explicitly specified you must use the names given. However, where not specified you are free to choose you own clear names for variables, functions, classes, files, and so on. The same applies to types and function signatures.

Your code must be able to be built from the top level folder (COMP0210Assignment1) using the following commands:

cmake -B build

cmake --build build

You can assume that when we build your project we will have the necessary dependencies

(FFTW3, Catch2) installed. It is important that you do not add further dependencies.

Installing FFTW3

If you are not using Docker you will need to install FFTW3 on your machine. To do so:

1. Download FFTW 3.3.10 from the FFTW website.

• You can also use the terminal command: wget http://fftw.org/fftw-3.3.10.tar.gz

• Unzip it with: tar -xvzf fftw-3.3.10.tar.gz

• Move into the folder: cd fftw-3.3.10

• Build and install:

– mkdir build

– cd build

– cmake ..

– make -j

– make install

You must follow these instructions to build and install using CMake, otherwise CMake will not be able to find FFTW3 when you try to build your project.

Reporting Errors

Contact the lecturer directly if you believe that there is a problem with the starting code, or the build process in general. To make error reporting useful, you should include in your description of the error: error messages, operating system version, compiler version, and an exact sequence of steps to reproduce the error. Questions regarding the clarity of the instructions are allowed. Any corrections or clarifications made after the initial release of this coursework will be written in the errata section and announced.

Assignment: Image Reconstruction with Optimisation Algorithms

In this assignment you will write a library for solving optimisation problems of a given form. This library will be designed to be highly generic, and make use of the design features that we have discussed in C++ and OOP. You will then write application code which uses this library to solve image reconstruction problems from simulations of lossy image data.

0.1 Mathematical Background

0.1.1 The Optimisation Problem

Optimisation problems generally involve finding inputs which minimise or maximise a given function, called a “cost function”, of those inputs. Many simple optimisation algorithms involve using the gradient of the function at a given point as a way of generating a new guess for a point closer to the minimum. In our case, we want to minimise a cost function of this form:

C(x) = f(x, y, Φ) + g(x)

where:

· x is the input that we want to vary in order to minimise C. In our case, this will be an image model, which we are trying to get as close as possible to the “true” image.

· y is a set of measurements. These are typically noisy, incomplete, and/or distorted in some way. For example, an image may be blurry or have data missing, or suffer from other camera effects that distort the image that it is trying to capture. Our goal is to reverse these processes as well as possible.

· Φ is a linear operator which translates an image model (x) into the measurements (y) that would be made from that image: y = Φ(x). (In this assignment we will ignore any random measurement noise, although in practice this is usually present!) Φ is sometimes known as a measurement operator, and can capture things like incompleteness of the image (patchesmissing) or camera effects such as point  spread. For a given linear operator Φ, we will also need its adjoint Φ † . The adjoint is not the same as the inverse, but it is often used as a pseudo-inverse (approximate inverse) for operators which are not invertible. Most linear operators that we are interested in are not invertible.

· f is a differentiable function, which is to say that the gradient ∇f is defined. A commonly use differentiable function in cost functions is a Gaussian log-likelihood L ∼ |Φ(x)−y| 2/2σ 2 . (You will not need to be able to differentiate the functions yourselves; we will provide the mathematical functions and their gradients.)

· g is a non-differentiable function, so ∇g is not defined, and hence we cannot rely on the gradient alone. Instead we can defined other ways of taking a step motion such as a proximal operator. In our case we will use this function to define a non-smooth prior, commonly found in optimisation problems.1

N.B. You do not need to understand the mathematics of gradients, proximals, or linear operators to implement the code, all you need to know is that differentiable and non-differentiable functions are treated separately with one defining a gradient and the other defining a proximal. The functions, gradients, and proximals will be given to you. You should however bear in mind what the input and output types of each of these kinds of functions are.

This optimisation problem is very general, and not limited to the application to imaging problems that we will look at, nor do f and g need to be a likelihood and prior respectively. To implement this generality we will use abstract and generic code, with increasingly specialised implementations as we work towards solving specific imaging problems. 

0.1.2 An Example Problem

Let’s consider that the real image that we are trying to recover is this:

The perfect image is however always unknown to us in real situations. Generally our measurements have imperfect information; let’s assume in this case that some of the data is missing. In this case the linear operator Φ would take an image x and return a subset of those pixels. This linear operator loses information, and so can’t be reversed perfectly without knowing the originalimage. This means that  the adjoint of our linear operator cannot reconstruct the image, and instead just places zeroes (black) in the missing pixels. The image that we start with from our measurements is therefore something like this:

From our own intuitive understanding of images we know that this these black pixels everywhere are very unlikely to be correct: these pixels should probably have similar values to pixels that neighbour them. Therefore, using a prior which encourages some kind of smoothness, we can reconstruct a better approximation of the image like the one below: 

This image has been reconstructed from the sub-sampled image using the algorithm described in this assignment. The cost function used to reconstruct this image needs to balance the smoothness prior with the fidelity to the data (the pixels that we know are correct shouldn’t change too much) to avoid over-smoothing.

0.1.3 Solving the Optimisation Problem: Iteration and Convergence

The solution to the optimisation problem will be a simple iterative algorithm. (This can be made faster and more robust by a number of improvements which are beyond the scope of this assignment!)

The algorithm will stop once either:

• a convergence criterion is met, or

• a maximum number of iterations is reached (to prevent infinite loops).

In each iteration the algorithm will take the estimate xn and produce a new estimate xn+1.

These are related by the formulae:

x 0 = xn − α∇f(xn)

xn+1 = x 0 − βproxg (x 0 )

where

· x 0 is a temporary intermediate value calculated by taking the gradient step

· proxg is the proximal operator of the function g

· α and β are constants which control how big the step taken is relative to the gradient / proximal step. Tweaking these parameters can have an impact on how dominant either function (f or g) is in determining the solution.

The convergence criteria are used to determine when to stop the iterative algorithm. We will explore two common convergence criteria:

This means that we are checking that the magnitude of fractional change in either x itself or the cost function C(x) in a given iteration is smaller than some tolerance δ.

These criteria will be expressed in the form of a function H(xn, xn−1, δ) which returns True if the convergence criteria are met (stop the algorithm) and False otherwise.

Start by familiarising yourself with the starter repository.

There are multiple folders which contain source files:

1. app is for application code. These are linked to the OptimisationLib.

2. Optimisation contains library code such as I/O, algorithms, functions / linear operators etc.

3. test for test code. There are subfolders for data and images that will be used during testing.

4. data contains data that you will use when running your final applications.

5. images is where you will place the image files from your reconstruction process.

You will also have a README.md and Responses.md file. The README.md should be filled in as described in the final section of this assignment, and should contain all the necessary information for a reader to build and run your application. Responses.md is used for answering written questions that appear in this assignment, such as discussing your design choices.

You will find in the source files:

1. Some partially coded tests that you will need to adapt to use your classes / data structures.

2. Methods for reading and writing data files (.dat files in the data folders) and image files (in .pgm format).

3. Additional utilities such as methods for std::vector arithmetic.

4. Some optimisation specific utility code that you will make use of throughout the assignment such as generating Gaussian kernels and a simple convergence function.

5. Some simplified example classes with the structure of the differentiable functions, non differentiable functions, and linear operators.

You should inspect the CMakeLists.txt files as well and understand the structure of the build system. This is necessary in order to understand how to add new files to the code base.

0.4 Writing In Responses.md

Please be sure to clearly label each question response with the section of the assignment that you are answering. Your file should be well formatted so that it is easy to read, you should write clearly, and explain your ideas as fully as you can.

1 Library Code Part 1: I/O [8 Marks]

This project will require file i/o for two main types of files: data files and image files. Data files will be used to load measurement sets (y), as well as some information related to linear operators. Image files are used specifically for images, although the format is very similar.

The format for data files is as follows:

• The first line contains the dimensions of the image to be reconstructed separated by a space.

• The second line contains the number of data in the set.

• The third line onwards contains the data, with one item per line. A data file may use any arbitrary type for its data, but all data in a given file must be the same type. In this assignment we will consider both real (floating point or integral types) and complex (std::complex<double>) types.

– A real value is just a single number.

– A complex value is represented by two numbers separated by a space. The first value is the real part, and the second is the imaginary part. (The imaginary part is strictly numeric, there is no i or j character to denote the imaginary part.)

Image files being read are in .pgm format. To view them you can install a viewing extension for VSCode or use an image program such as gimp.

We have provided functions to do the following:

• Read and write well formatted data files for real and complex data.

• Read and write image files in .pgm format.

Tasks

1. Add error handling to the data file reader.

2. Test your data reader function.

N.B. You do not need to also add error handling to image reader, or to the write functions, but you may if you wish. (The error cases for images and data are very similar.) The image read and write functions are limited compared to the flexibility of the .pgm format, and only designed to read and write with specific layouts in mind, so you may need to modify if you wish to use them with images that we have not provided. [4 marks implementation; 4 marks testing]

2 Library Code Part 2: The Generic Problem [12 marks]

2.1 Abstract Classes [4 marks]

Recall the definition of the cost function, which we are trying to minimise:

C(x, Φ, y) = f(x, Φ, y) + g(x).

The three main components that define a given optimisation problem are the differentiable function f, the non-differentiable function g, and the linear operator Φ. In order to describe theproblem in abstract terms, we don’t need to know exactly what any of these functions actuallyare, we only need to know how they can be used. This makes them a good candidate for abstractclasses.

In this section you will write an abstract class for each of these three types of function. You may use the Quadratic, Identity, and Indicator classes as guidance for the kind of functionality that they need to perform, but bear and in mind that these are simple, concrete cases, and you should be writing classes which are abstract and generic.

2.1.1 Differentiable Functions

A representation of a differentiable function f needs to be able to do two things:

• Evaluate the function for a given input, f(x). This output is a single real number (double).

• Evaluate the gradient for a given input, ∇f(x). This output is a vector the same size and type as x.

Tasks

1. Write an abstract class to implement a differentiable function.

• The type of the data in the input and the gradient should be generic, and must be the same.

• The return type of the function evaluation should be a single double.

• The function evaluation should overload operator().

• The gradient should be a function called gradient.

2.1.2 Non-Differentiable Functions

A representation of a non-differentiable function g needs to be able to do two things:

• Evaluate the function for a given input, g(x). This output is a single real number.

• Evaluate the proximal operator for a given input, proxg (x). This output is a vector the same size and type as x.

Tasks

1. Write an abstract class to implement a non-differentiable function.

• The type of the data in the input and the gradient should be generic, and must be the same.

• The return type of the function evaluation should be a single double.

• The function evaluation should overload operator().

• The proximal operator should be a function called proximal

2.1.3 Linear Operators

To implement linear operators in our code we need to represent both the linear operator itself (sometimes referred to as the forward transform for clarity) and its adjoint.

· The linear operator acts on input vectors x and produces measurement vectors y. The types of data contained in x and y do not have to be the same. For example a linear operator might take a vector of real numbers and return a vector of complex numbers.

· The adjoint operator acts on measurement vectors y and produces input vectors x. This is used for obtaining the initial guess, amongst other things.

Tasks

1. Write an abstract class to represent a linear operator.

• The type of data in the input and measurement vectors should be generic.

• The forward transform, Φ(x), should be implemented by overloading operator().

• The adjoint , Φ † (y) should be a function called adjoint.

2.2 The Iterative Algorithm [5 marks]

The details of the algorithm are covered in the mathematical background section 0.1.3 of this assignment.

Tasks

You should now implement the iterative optimisation algorithm as a function.

1. This function should return the final estimate of x after meeting the convergence criterion.

2. The function should take the following parameters:

• f(x), the differentiable function

• g(x), the non-differentiable function

• y, the measurement set

• Φ, the linear operator

• α, the step size for the gradient

• β, the step size of the proximal

• nmax, the maximum number of iterations before giving up

• δ, the convergence threshold

• H(xn, xn+1, δ), the convergence function (optional with default; see below)

3. The algorithm should terminate if H(xn, xn+1, δ) returns True or if the maximum number of iterations is reached.

4. The function should generate an initial guess by using the formula x0 = Φ†y before the first iteration.

• The result of the first iteration is therefore x1, and if convergence is not reached then the result returned will be xnmax .

5. Your algorithm function should be programmed to work with the generic interfaces that you have already developed in the sections above.

6. It should be able to accept an arbitrary convergence function H as long as it matches the correct signature, as indicated by the provided convergence function bool converged(vector<Tx> &x1, vector<Tx> &x2, double tol) (see OptimisationUtils.h). This convergence function should be the default if no convergence function is given as an argument.

7. When deciding the types for this algorithm you should aim to keep it as generic as possible while respecting type constraints of the problem.

2.3 A Simple Test [3 marks]

Tasks

In this test we will check that we can solve a very simple problem: find the value x which minimises the function x 2 . In this case:

• We can use f(x) = x 2 (using the Quadratic class).

• We can use Φ(x) = x (using the Identity class).

• We can set g(x) = 0 to make no contribution (using the Empty class).

In this case:

x0 = Φ† (y) = y, so the value of y is your initial guess.

Each iteration of the algorithm will take a gradient step ∆x = 2αx.

1. Modify the Quadratic, Identity, and Empty classes so that they are usable with your optimisation algorithm.

2. Use the default convergence criterion.

3. The test will check that your result differs from 0 by no more than 10−3 .

• How close to 0 you get will depend heavily on your convergence tolerance, so you may need to experiment to set your tolerance low enough.

• α should be small enough that you don’t immediately over-shoot your target. How you set α will affect the number of iterations your test needs.

4. In Responses.md describe some limitations of this as a test of our algorithm.

3 Library Code Part 2: Concrete Linear Operators [21 marks]

In order to try our code out on a meaningful imaging problem, we will create some more interesting linear operators and functions. In this section you’ll implement two linear operators:

• A subsampler, which takes a sub-set of data from an image to produce a new vector with less information.

• A convolution operator, which performs a convolution. A common application of convo lution in imaging is point spread which models how a single point in the real world gets rendered by a system such as a camera. One example of this would be Gaussian blur, a common model for blurring images.

3.1 Implementing a Sub-Sampling Linear Operator [7 marks]

A sub-sampling operator receives a vector of data of a given size, and returns a sub-set of that data. This linear operator can simulate data loss, or patchy data (for example incomplete detector coverage of an area).

In the following we will use the notation I to represent the list of indices preserved by the subsampling process. A sub-sampling operator is determined by the size of image it is intended to work on, and this list of indices.

Tasks

Write a sub-sampler linear operator class such that:

1. The list of preserved indices I should be passed into the constructor.

2. The forward transform takes an input x and returns a measurement y containing all the input values at the indices given in the constructor. For example, if the linear operator is constructed with the list of indices I = [1, 4, 2], then when applied to a vector x = [0, 9, 4, 1, 2, 5] it should return the smaller set y = Φ(x) = [9, 2, 4].

· The list of indices does not need to be ordered, and the output vector should contain the elements of the input vector in the order given in the list of indices I. It isimportant that you respect the ordering in I because some problems that you will runwill involve reconstructing an image from a sub-sampled image and a list of indices.

3. The adjoint takes a measurement set y and returns an image x 0 with:

· If y = Φ(x) then the values in x 0 = Φ† (y) should match those in x for all preserved indices i.e. xi 0 = xi for all i ∈ I.

· The values at all other indices, where data has been lost, should be filled with 0.

· Using the same example above, where I = [1, 4, 2], x = [0, 9, 4, 1, 2, 5], and y = Φ(x) = [9, 2, 4], then x 0 = Φ† (y) = [0, 9, 4, 0, 2, 0].

4. You should consider error handling carefully when implementing this class.

5. Write unit tests for your sub-sampling linear operator.

3.2 Implementing a Convolution Linear Operator [14 marks]

N.B. In this section you will need to use the FFTW library to implement convolutions using Fourier transforms. For more information on this library and its syntax please check the FFTW section (0.2) at the end of the introduction.

In this example you will create a more complex linear operator. A convolution is a common operation on images which can simulate things like blurring and lens effects. Although convolutions can be calculated directly, we will calculate it in a slightly different way using a Fourier transform, as this makes it easier to define the adjoint operator. You don’t need to understand exactly what a Fourier transform (or convolution) does, except for the following:

• A Fourier transform converts between two representations of a function. An image x can be considered a real space function, where the function is given a value at various points in space. A Fourier transform gives the Fourier space function, which we notate with a tilde: F T(x) = ˜x. The Fourier space function represent a function as a sum of waves, with amplitudes defined at various frequencies instead of points in space. (For example, a Fourier transform would convert a sin wave in real space to a single peak in Fourier space at the frequency of that wave.)

• The Fourier transform has an inverse, which converts from Fourier space to real space: F T −1 (˜x) = x.

• Given an image x and a kernel k, we can calculate the convolution x 0 = x ∗ k by the following process:

– Apply the Fourier transform to obtain ˜x and k˜.

– Calculate each element of ˜x 0 by a pointwise multiplication: ˜xi 0 = ˜xik˜i .

– Apply the inverse Fourier transform to get x 0 .

The benefit of this approach is that it makes the adjoint (the “deconvolution”, which is typically extremely difficult to calculate in real space) fairly simple to define.

• Given a convolved image x 0 and a kernel k, we can (approximately) reconstruct the original image x by the following process:

– Apply the Fourier transform to obtain ˜x 0 and k˜.

– Calculate each element of ˜x by a pointwise division: ˜xi =

– The deconvolution process can often be unstable, especially if some of the kernel values k˜i are small. It is more stable to implement the following modification: ˜xi =  for a small  to prevent the denominator from being too small.

– Apply the inverse Fourier transform to get x 0 .

N.B. Because the convolution is implemented as a point-wise multiplication in Fourier space, the dimensions of the image buffer and the kernel buffer must be the same.

Tasks

Using the definitions above, write a convolution linear operator class which implements the following:

1. The constructor should take an image representing the kernel. The kernel must be the same size as any images which get passed to this operator.

2. The forward transform should perform a convolution. It should take a vector of real data representing an image and return another vector of real data.

3. The adjoint should be the deconvolution. It should likewise take an image and return an image.

• You should use  = 0.01 by default in the pointwise division to match the exampledata given for testing.

• You should take  as an optional argument in the constructor.

4. The image input and output data must be real (double), but the intermediate Fourier space representations ˜x 0 and k˜0 must be complex (fftw_complex*).

• When applying the forward transform to an image you should take the imaginarypart of every pixel to be 0.

• After applying the backward transform you can discard the imaginary part to get areal image.

5. You should use the FFTW library to perform the forward and backward fourier transforms.

6. Test your convolution class by using the test images provided.

• The outline of a test for this operator has been provided for you.

• Bear in mind when comparing your images to the ones provided that the values in the image files are rounded to the nearest integer.

Hints

• An object of this class should be able to be applied to any supplied image of the correct dimensions, i.e. the source image should not be fixed at construction.

• You should try not to do unnecessary work inside the class as this operator will be calledevery iteration in your optimisation algorithm.

• Don’t forget the scaling on the inverse transform!

4 Library Code Part 4: Concrete Functions [22 marks]

4.1 Implementing a Gaussian Likelihood function [6 marks]

In this section you will implement an non-normalised Gaussian log-likelihood function as anexample of a differentiable function. The functional form is: f(x, y, Φ) = (2σ 2 ) −1 |Φ(x) − y| 2

• This likelihood is a measure of the difference between the measurement data that we have (y) and the measurements that we would expect from our image after applying themeasurement operator (Φ(x)).

• The parameter σ is the standard deviation of the likelihood, and it determines how sensitiveour function is to variations in Φ(x).

• In our imaging applications the image, x, is always a vector of real values.

• The measurements, y, may be a vector of real values or a vector of complex values,depending on the measurement operator.

– For example, if Φ is a simple sub-sampler, then Φ(x) and y are vectors of real values. But if Φ were a Fourier transform, then Φ(x) and y would be complex, because the output of a Fourier transform is complex.

Note that the calculation for the magnitude of a vector |v| 2 is different for real numbers and complex numbers. For real numbers we have:

|v| 2 = P vi 2 ,

and for complex numbers we have:

|v| 2 = P vi × vi ∗ ,

where vi ∗ denotes complex conjugation. For a complex number (a + ib), its complex conjugate is defined as:

(a + ib) ∗ = (a − ib)

N.B. The C++ standard library provides the complex type std::complex<double>, which has the member functions real(), imag(), and conj() for getting the real part, imaginary part, and complex conjugate of a complex number respectively. The operators +, -, *, and / are also overloaded for complex types.

The gradient of the log-likelihood can be calculated as follows:

 

Tasks

1. Implement the Gaussian likelihood as a sub-class of the differentiable function class.

2. The constructor should take the measurement operator, the measurements, and σ as arguments.

3. Your class should be templated to be able to use real (double) or complex (std::complex<double>) types for the measurements, depending on the type of measurements the measurement operator will produce.

4. You should test your log-likelihood function and gradient following the guidance in thetest file.

4.2 Implementing an DCT L1-Norm Prior [9 marks]

Our non-differentiable function will make use of the discrete cosine transform (DCT), also provided by the FFTW library. (For more information, please see the FFTW section of the introduction.) Unlike the FT, the input and output buffers for this transform are both real. This transform is very common in image processing, such as jpeg compression. The DCT converts an image from a collection of pixels with different intensities to a collection of cosine-wave patterns with different intensities. The sum of all these cosine patterns gives you the image again.

For our purposes, the L1-Norm is simply a name for the sum of absolute values of the elements in the vector. By adding it to our cost function, we are penalising our solution for having large values in the DCT representation of the image. This kind of condition is often applied to encourage smoothness in an image.

Calculating the DCT L1-Norm involves three steps:

• Apply the DCT to the input image x.

• Sum the absolute value of each element in DCT(x).

• Return the sum.

Instead of a gradient, we have a proximal operator. In the case of the L1-Norm, the proximal operator is a soft thresholding. To calculate the proximal operator one must:

• Apply DCT to input image x.

• For each element in DCT(x), apply the soft thresholding function S with threshold τ :

– if ˜xi < τ then S(˜xi) = 0

– otherwise S(˜xi) = ˜xi(1 − |x˜ τ

• After applying the soft thresholding to the frequency space representation ˜x, apply the inverse DCT to get the modified image.

• Rescale modified image to account for unnormalised DCT.

• Return the difference between the original image and the modified image.

Tasks

1. Implement a DCT L1-Norm class as a non-differentiable function.

2. The constructor should optionally take the threshold parameter τ , which should default to 0.1. In order to set up the DCT, the constructor will also need the dimensions of the image.

3. Use the FFTW library to perform the DCT and its inverse.

4. Implement both the function and the proximal operator.

5. Test your class using the test framework provided.

• N.B. You should use the updated TestSphere_proximal.dat provided if you have used the corrected FFTW interface from this version of the assignment text.

Hints

• As with your convolution class, an object of this class should be able to be applied to any supplied image of the correct dimensions, i.e. the source image should not be fixed at construction.

• You should try not to do unnecessary work inside the class as this operator will be called every iteration in your optimisation algorithm.

4.3 Convergence Functions [3 marks]

The last part of our algorithm that we need in order to run an iteration problem is a convergence function.

We have provided one convergence function in OptimisationUtils.h, which looks like this:

H(xn, xn+1, δ) = |xn+1 |x −xn|

n

| < δ.

Tasks

1. You should implement the following convergence function:

H(xn, xn+1, f, g, δ) = |C(xn C +1 ( ) x − ) C(x)| < δ.

2. Explain in Responses.md how you can ensure this function has the correct signature to be passed into your optimisation algorithm once you have defined all the pieces of your optimisation problem. Give an explicit code example of how you would do this.

N.B. There are no marks for testing this function.

4.4 Integration Test [4 marks]

You should now perform a check that your overall algorithm is working with all the componentsin place. You should use the image test/data/images/TestSphere.pgm as the basis for yourtest; this will be the “true” image that you try to recover. You can simulate the measurementprocess by applying a measurement operator to the true image to generate a set of measurements.

Tasks

1. Write a test which checks the performance of your algorithm with your Gaussian Likelihood function, a DCT L1-Norm, and a sub-sampler which discards 20% of the data.

• You can use the following parameters:

– set σ = 0.1, the standard deviation of the log-likelihood.

– set α = 0.01, the gradient step factor.

– set β = 104 , the proximal step factor.

– set δ = 0.4, the convergence tolerance. (If you function converges too quickly,lower the tolerance. Don’t worry too much if it doesn’t converge, as long as itreturns a good result at max iterations. You may need to run for up to 50-100iterations.)

– set τ = 0.1, the soft threshold in the DCT L1-Norm.

2. You should check that reconstructed image that you get is an improvement on the dirtyimage (x0) by comparing both to the true image.

5 App Code Part 1: Simple Imaging Problems [10 marks]

In this section you will use the concrete classes you have created so far to solve two imagingproblems. We have provided the measurements data for each of these problems, but not theoriginal images. You must write an application which can reconstruct these two images, andsave your results to the images folder.

• One data set is an image which has been sub-sampled.

• The other is an image which has beeen blurred using a convolution with a sinc kernel.

Both problems will use a Gaussian log likelihood for f and a DCT L1-Norm for g, but they require different linear operators.

Tasks

1. Write a command line application called ImageReconstructor. You should implement the following command line arguments:

• -h: if this flag is used then the application should display a help message and then terminate; the optimisation algorithm should not be run. If additional options are specified then they should be ignored.

• -f: the measurements data file

• -alpha: the value of α, the gradient step size

• -beta: the value of β, the proximal step size

• -sigma: the value of σ in the log-likelihood function

• -delta: the value of δ in the convergence function

• -k: the Gaussian kernel size in pixel (one standard deviation) only for a convolution problem

• -i: the indices file only for a subsampling problem

– The -k and -i flags should be mutually exclusive, but the user must supply one.

2. If the kernel option is provided then the app should create a convolution linear operator using a Gaussian kernel of the given standard deviation.

• You can generate a kernel of a given image size using GenSincKernel from ImageUtils.h.

3. If an indices file is provided then the app should create a sub-sampling linear operator using the provided list of indices.

• The indices can be loaded using the ReadData function in OptimisationUtils.h.

4. If the arguments are valid then the program should do the following:

• Calculate the dirty image (Φ † (y)) and save it to a file. This should be named the same as the measurement input file name with _dirty appended before the file extension. (The file extension should be .pgm.)

• Run the optimisation problem using the appropriate differentiable function, non differentiable function, linear operator, and parameters.

– σ = 1.0

– α = 1.0

– β = 104

– τ = 0.1

– max steps = 500

• Save the reconstructed image to a file. This should be named the same as the measurement input file name with _reconstructed appended before the file extension.

(The file extension should be .pgm.)

• Your solution may not converge before reaching the maximum number of iterations.

This is not necessarily a problem and your optimisation algorithm should still return a result which you can use.

5. Report all the command line parameters used for each run in your Responses.md file. Run your app for the following two problems:

• A sub sampling problem using:

– the measurements file data/measurements/UtahTeapot_subSampledMeasurements.dat

– the indices file data/operators/UtahTeapot_sampledPixels.dat.

• A convolution problem using:

– the measurements file data/measurements/UtahTeapot_convolved.dat

– a sinc kernel with scale 1.0.

6. Make sure that you place your dirty and reconstructed images into the images folder before you submit.

Hints

• You can use basic C++ command line parsing: http://www.cplusplus.com/articles/ DEN36Up4/ or you could try integrating a command line parser such as: https://github. com/CLIUtils/CLI11 (which is header only). If you want to use CLI11 you should include it as a sub-folder inside your submission and make sure that you modify your cmake so that your project builds correctly as the markers will not necessarily have this library on their own machine.

• You should compile with optimisation flags -O2 or -O3 turned on in your cmake and debugging turned off (no -g flag) to speed up running these problems.

• The quality of reconstruction can be sensitive to parameters, so you will not be marked on the quality of your final images as long as your code is functioning properly.

• The sub-sampling problem should show a clear improvement. The blurring problem may have a final result similar to the dirty image depending on how small the value of  is and how accurate the deconvolution is. You should see a clear improvement if  = 0.1.

• The parameters suggested will not necessarily give the best results, but should work without the solution diverging.

6 Library Code Part 4: Further Operators and Abstract Operator Composition [12 marks]

In this final part of our library development, we will make our linear operator code more powerful allowing new linear operators to be defined as the composition of two linear operators. We will also introduce one final linear operator, the Fourier transform operator.

6.1 Fourier Transform Operator [3 marks]