Quantum median filter for Total Variation image denoising

In this new computing paradigm, named quantum computing, researchers from all over the world are taking their first steps in designing quantum circuits for image processing, through a difficult process of knowledge transfer. This effort is named Quantum Image Processing, an emerging research field pushed by powerful parallel computing capabilities of quantum computers. This work goes in this direction and proposes the challenging development of a powerful method of image denoising, such as the Total Variation (TV) model, in a quantum environment. The proposed Quantum TV is described and its sub-components are analysed. Despite the natural limitations of the current capabilities of quantum devices, the experimental results show a competitive denoising performance compared to the classical variational TV counterpart.


Introduction
Digital image processing is an extremely, computationally demanding, strategic research field. This has led, over the years, to the development of many complex image processing algorithms on highly parallel, specialized hardware platforms. With the rapid progress of parallel hardware, suitably high performance is now available also for sophisticated imaging tasks. In the present study we focus on image denoising, an essential requirement for any other image processing, which refers to the recovery of a clean sharp image from a noisy observation.
In [21] Rudin, Osher and Fatemi presented one of the main mathematical models and algorithms for image denoising, the so-called Total Variation denoising (also known as TV model ). Since then, TV model has become a quite popular technique, which usage allows to improve overall image quality when the images are affected by noise or corruption, while well preserving edges and details. Efficient solutions have been proposed over time for the numerical optimizations of the TV model [3,4], but they are not able to fully exploit parallel computational resources. In [13], the authors designed a new optimization algorithm which is simple and highly parallelizable, and relies on median value computations, thus reducing computational effort to a sorting problem. Thanks to its low complexity, this algorithm is prone to be implemented on low-end devices or, more generally, in situations where a reduced amount of resources is available. That is the case of Quantum Image Processing (QIP), a novel and promising research field which goal is the development of image processing techniques for quantum computers, exploiting peculiar features from quantum world, like entanglement, superposition, interference [25]. Quantum computing is universally acknowledged for its ability to process data providing computational speedups compared to the classical paradigm. However, image processing is actually one of the most demanding applications in terms of resources for quantum computers. As a consequence QIP is still in its early stage, and thus is facing several fundamental problems, such as how to represent and store an image in quantum computers appropriately, and how to efficiently implement image processing algorithms [10,16].
In the present days, quantum devices are subject to several issues (i.e. noise, absence of error correction, low amount of available qubits, etc.), which defines what researchers commonly refer to as Noise intermediate-scale quantum era, or just NISQ [19]. This is an intermediate development phase where, through the exhibition of quantum devices' limits, algorithms have been designed to be as optimized and simple as possible [2]. In this context, QIP is a complex subject to deal with because several difficulties may arise in the attempt of implementing the classical image processing algorithms on quantum devices (some of them will be discussed later in this work). Many proposals have been submitted during the past few years [23], in the attempt to provide QIP algorithms that fulfill NISQ requirements. QIP algorithms are usually tested using quantum simulators or similar execution environments.
Here along this line, the presented work aims to solve the TV model in a quantum environment. The result is a Quantum TV filter, which integrates a Quantum Median Filter proposal by Li et al. [10,14]. This work focuses not only on the development of a QIP technique for image denoising which mimics its variational TV counterpart, but also on presenting an useful review of detected problems in the newly developed QIP field, and offering possible solutions and ideas for future developments and optimizations.
The paper is organized as follows. In Section 2 we introduce basic theory about Total Variation denoising technique and a median formula for efficiently solving an anisotropic TV problem. In Section 3 we discuss Quantum Image Processing concepts and related issues. In Section 4 we explain in detail how Quantum TV Filter works, highlighting its quantum components, and analysing its circuit complexity in Section 5. In Section 6 we focus on the analysis of experimental results, comparing both quantum and classical implementations of the TV algorithm. Section 7 reports conclusions and future works. In the Appendix we briefly report details on a few quantum modules used in the design of the proposed Quantum TV.
For basic notations and fundamental knowledge about quantum image processing computing, we refer the reader to [25].

Total Variation image denoising
The goal of denoising is to obtain an image u * not only with small variations in intensity between pixels but also close to the observation f . At this aim, the class of variational methods for image restoration relies on determining restored images u * ∈ R N , given a noisy image f ∈ R N , as the minimizers of suitable cost functionals J : R N → R such that, typically, restoration is casted as an optimization problem of the form: where the functionals R(u) and F (u; f ), commonly referred to as the regularization and the fidelity term, encode prior information on the clean image u and the observation model, respectively, with the so-called regularization parameter λ > 0 controlling the trade-off between the two terms. In particular, the functional form of the fidelity term is strictly connected to the characteristics of the noise corruption. A classical choice for the fidelity measures the data fitting in terms of the 2 -norm, in formulas: For what regards the regularization term R(u) in (1), a very popular choice is represented by the Total Variation, presented in the following two forms: where D x , D y denote the horizontal and vertical operators, respectively, and D = (D x , D y ) is the gradient operator ∇ in the discrete setting. By substituting the TV regularizer T V (u) := R iso (u) in (3) or T V (u) := R ani (u) in (4) and the fidelity term (2) for R and F in (1), respectively, one obtains the so-called TV-L 2 -restoration model, originally introduced in [21]. In formulas: (5)

Median formula for TV
Li and Osher in [13] proposed an efficient and highly parallelizable method for solving TV model (5) with anisotropic TV regularization (4). They proposed to solve iteratively N one-dimensional optimization problems to obtain an accurate solution of the N -dimensional optimization problem (5). In particular, for each pixel u ∈ R, we consider the local minimization problem: where F (u) = λ(f − u) 2 , u * , f ∈ R are respectively the denoised and noisy versions of the same pixel, while u i belongs to the set of κ neighboring pixels and w i ≥ 0 are given weights. In [13], a simple method for computing u * in (6) is derived and here reported for self-consistency.
Theorem 1 Supposing the w i are non-negative and the u i are sorted as u 1 ≤ u 2 ≤ ... ≤ uκ, the function F is strictly convex and differentiable and F is bijective; then the minimizer of (6) is a median: where p i = (F ) −1 (W i ) and In our formulation, the neighborhood of the current pixel u are simply u u , u d , u l , u r , the vertical and horizontal direct neighbors pixels, respectively. The adopted 4-neighbors strategy allows us to apply the median formula in parallel on multiple pixels at one time using a proper configuration, since each pixel is directly affected only by its 4 neighbors.
The fidelity is defined as F (u) = λ(f − u) 2 , and consequently where W is a sum of weights previously defined in (8). The denoised pixel u * is then obtained as the median value: where the p-values p i are calculated following Theorem 1, with w i = 1 and κ = 4: The image denoising problem can hence be solved by iteratively computing (7) pixel-by-pixel over the whole image until convergence, which is guaranteed by the following result.
Theorem 2 The algorithm defined by repeatedly applying (7) at the jth pixel, converges, i.e. u This means that, after k iterations, with k → ∞, we obtain the minimizer of problem (5).
For the numerical implementation the stopping criterion considered is computed as follows: given a small tolerance value , process stops when, at iteration k where u (k−1) and u (k) are consecutive processed images. The resulting algorithm is described in Algorithm 1.

Algorithm 1 Anisotropic TV algorithm
3 Quantum Image Processing QIP aims to design quantum algorithms that, once constructed a quantum state which encodes an image, implement image processing techniques in a quantum environment. QIP is a novel topic in quantum computing and many issues are far from being solved, as explained in [20]. A central issue regards the image encoding in a quantum environment, which will be discussed in Section 3.1. Another limiting issue is that, nowadays, QIP represents a demanding branch of quantum computing, as it needs a lot of resources that are far from being offered by current or medium term devices. At the moment, researchers are proposing many different approaches to the problem, even though only a limited number of these are commonly accepted and used [23][24][25]. In order to tackle the memory restrictions, we subdivide the image into a set of sub-images as described in Section 3.2.
The denoising QIP here proposed consists of three steps, illustrated in Fig.1. The image is first encoded in the quantum environment, then it is processed by quantum circuits to perform the denoising task, and finally the denoised image is measured to convert it in a classical image format. The used measurement process is detailed in Section 3.3.

Quantum Image Representation
A quantum image encoding is defined Quantum Image Representation (QIR). Unlike classical image processing, where a set of well-known standard formats are available and well-assessed, in QIP many encoding QIR techniques were proposed and tested [14], but on the other hand none nowadays has distinguished itself as standard. The most used QIR technique is the Novel Enhanced Quantum Representation (NEQR).
This kind of representation needs to encode the following image's data: • Pixel coordinates, encoded by qubits |XY . An image of dimension D X × D Y , usually needs to use d X = log 2 D X qubits for horizontal coordinates and d Y = log 2 D Y qubits for vertical ones. In this way • Pixel value, encoded by one or more qubits |C . It uses q qubits for representing N q = 2 q possible values in a binary encoding.
Without loss of generality, we consider grayscale images I with a square 2 n × 2 n domain and values in the range [0, 2 q − 1]. In this way we consider |XY with 2n qubits, where n = log 2 D , and |C using q = log 2 N q qubits to represent 2 q possible values. Therefore, a grayscale image I is encoded in the following quantum state of 2n + q qubits: where and each C i XY qubit is |0 or |1 . A simple example of NEQR representation for a 2 × 2 grayscale image is illustrated in Fig. 2(left), where each corner number denotes pixel's coordinates, while centered value indicates pixel's intensity [25]. The quantum wave function for this image in NEQR is hence the following: This is obtained by putting the coordinate qubits |X and |Y in a superposition state using Hadamard gates, then entangle them with qubits |C XY using a series of CNOT gates. The resulting NEQR circuit is illustrated in Fig. 2(right) for the image in Fig. 2 NEQR is a very versatile representation for image computation. However, it presents some drawbacks. As pixels are encoded one by one, large images produce long NEQR circuits; this means that circuit depth grows linearly with image size. As the image size grows, the number of coordinate qubits become larger and the control qubits needed for encoding quickly become more than two. A Toffoli gate with more than two control qubits is defined Multiple CNOT gate, or just MCX: this operator is decomposed before execution, using many Toffoli gates and some auxiliary qubits called ancilla. From an efficiency point of view, more coordinate qubits means larger MCX to be used, which leads to a larger amount of gates.

Image pre-processing
In order to reduce the memory usage, we split a pre-padded image (onepixel bordered) into many smaller overlapped patches, 4 × 4 sub-images. Once extracted, a patch of pixels will be processed by the quantum TV algorithm. At the end of each iteration, will be necessary to reassemble the resulting image from the output patches. Fig. 3 illustrates the image pre-processing procedure for a sample image: red square frames original image, while green squares in patches highlight the processed pixels.
In order to speed up quantum circuits generation, we have subdivided workload using multi-threading: each thread is tasked to assemble a QTV for each image patch, using pre-assembled circuits and generating remaining ones. This approach considerably accelerated the execution process for what regards the generation phase.

Image measurement
The image extraction from the QIR format is a not-trivial process and its performance depends on the quantum representation used in the algorithm.
A quantum representation collects all image's data in a single quantum state. Due to the nature of this particular quantum state, the extraction is not a deterministic process: for each measurement, one of the possible pixel coordinate-value association is randomly obtained as outcome from the collapse of the quantum state. For a complete recovery of an image, it is necessary to execute the same algorithm many times.
The image measurement is the last step in Fig.1.

Quantum TV Filter
In this section we introduce the Quantum TV Filter, named QTV, a quantum implementation of the TV regularization algorithm described in Section 2 applied to qubits involved in the NEQR quantum image representation. This work extends and improves the work in [10] where the authors proposed a quantum solution for implementing a simple median filter for image processing.
The proposed circuit processes an image input in NEQR representation and provides a denoised image in output in the same NEQR form, according to the scheme in Fig. 4(left).
According to the TV algorithm presented in Section 2 we have to iterate a core process for each pixel of the input image. Considering a four-pixel neighborhood configuration, this process is composed by three steps defined as three different quantum operators acting for each pixel:  The quantum operators are assembled into the QTV structure illustrated in Fig. 4(right).
In the following, we will describe in detail the three operands which characterize the Quantum TV. Each one is composed of many sub-circuits, or modules.

Neighborhood Preparation
This operand is in charge of extract neighboring pixels from NEQR representation. At this aim, we used Cycle-Shift (CS) module to change a coordinate register's value, allowing us to shift an image up, down, left or right; see Appendix A for details.
The structure of the NP operand, which output is a quantum superposition state of the central f , up u u , down u d , left u l and right u r pixel values, is illustrated in Fig. 5. Specifically, once obtained a pixel value from NEQR, we use CS to shift coordinate values, then we re-apply NEQR to gather a new pixel value corresponding to the new position data. With the exception of the first NEQR for image encoding, we avoid the usage of H-gates applied to coordinate registers, as we want to extract a specific (and not random) color outcome.
The output of the NP operand is a quantum state obtained starting from |ψ = |0 ⊗(q+2n+5q) , which reads as |u l |u u |u r |u d |f |0 |X |Y (15) where each qubit in the |C register is reset to the |0 state. More in detail, following the scheme in Fig.5, we have

P-values computation
Starting from the obtained neighborhood values, we compute the p-values according to relations (10). This reduces to adding some constant values to f , the central pixel.
However the QTV algorithm only works in unsigned integer arithmetic, while the TV algorithm works with floating point numbers, thus getting more accurate and precise results. Therefore, instead of the p-values p i = f + Wi 2λ , we force approximated rounded values The p-values are then given by: The final design of the P-values Computation module is hence illustrated in Fig. 6, and it is composed by three sub-circuits for setting, adding and subtracting the mentioned constants. In particular: • SETTER module: assign the round to the nearest integer of a given value to a register, which is used to encode our constants; • ADDER module: add two values encoded in two quantum registers. We refer to the Appendix A for more details; • SUBTRACTOR module: subtract two values by using an adder module, according to α − β ⇒ α + β.

Median function
To determine a median from a set of values, we need to sort them. The sorting strategy here adopted follows the proposal in [10] for a set of 9 sortable elements. If we re-arrange these values in a matrix form, then we simply need to follow these three steps in pipeline: sort each column, sort each row, and sort the right diagonal. This will guarantee us to store the median value in the central quantum register. The strategy is illustrated in Fig. 7. The corresponding quantum circuit is shown in Fig. 8. The Sort module is then the core of Median Function module. Sort module has to order its three input registers. The total ordering of three elements is a trivial problem, as it is reduced to compare two positive integer values Fig. 7 The basic sorting strategy. and, if not ordered, swap them. A sequence of three comparisons is necessary and sufficient to reach a correct ordering. A Sort module is then a sequence of three sub-circuits, named Swapper (SWPR). A Swapper module is in turn composed by two sub-circuits: a Comparator (COMP) and a Controlled Swap (C-SWAP). Sort and Swapper modules are shown in Fig. 9.
The Comparator module evaluates two register values a, b and provides, on an auxiliary qubit e, the result of the comparison a > b, as follows This result will be next used to control a C-SWAP gate, so that if e = |1 , then the values a and b are swapped.
The Comparator module is described in Appendix A. The designed quantum circuit for the Comparator module is more efficient with respect to other proposals. For example, in case q = 8, this circuit uses less quantum elementary gates than other existing methods: depth = 64 with respect to the Sort module applied in [10] which has depth = 1.091.767.

Circuit complexity analysis
In order to estimate a quantum circuit efficiency, we have to look at its depth, that is the longest path in it. The path length is always an integer number, representing the number of gates it has to execute in that path. At this aim, we are going to analyze each operand of Quantum TV to derive an approximate estimation of its depth.
Neighborhood Preparation is for sure the most demanding operand of all filtering algorithm, because it uses multiple instances of NEQR circuit, which length is variable according to the number of pixel to encode. Considering a NEQR implementation that commits the least amount of MCX gates and encodes N pixels, its depth grows polynomially with N as follows, Other modules involved are Swap and Cycle Shift. Swap depth depends on the number q of color qubits, thus its depth is always equal to that value. Cycle Shift depth, instead, depends on coordinate register size and it is equal to log 2 √ N · (log 2 √ N − 1). Neighborhood preparation uses NEQR, SWAP and CS five times in a row. Hence the NP module overall depth is polynomial in N and we can estimate the NP depth as: For what concerns the P-values Computation depth, the setting phase involves four SET modules, which however are applied simultaneously, so they count as a single one. An Adder module depth instead is dependent on the q value. Full-Adder is composed of q Half-Adder, with constant depth (considering also Reset gates); then it is followed by 4q + 1 sequentially placed gates. Thus we have: ADD depth = 9q + 4q + 1 = 15q + 1 SUB depth = ADD depth + 2 = 15q + 3 From these results, we can derive the P-values Computation module total depth which is polynomial in q: To compute Median Filter depth, we just need to estimate SWPR module depth and count its occurrences in the circuit. To do so, we add up Comparator and C-SWAP depths.
Comparator's depth depends on q, as for each color qubit it uses 6 gates. Although Toffoli gate used in this module counts as three, as two additional X-gates need to be added in order to work. This means that module depth is estimated as 8q. C-SWAP depth is the same as SWAP. This means that Swapper's total depth is estimated as 8q + q = 9q.
In Median Filter there are multiple occurrences of SWPR module, which can be further reduced if row and column sorting are executed simultaneously: From this analysis we have derived that Quantum TV algorithm implements the NP module with polynomial complexity in the number of pixels N . The PC and MF modules instead have a polynomial complexity in q = log 2 N q , i.e. the logarithm of the number of colors, hence providing an exponential speedup.

Experimental results
In this section we evaluate the performance of the quantum TV algorithm (QTV) on denoising grayscale images, and we present some preliminary results from the comparison with the variational anisotropic TV in Algorithm 1.
The reference images used for the test have dimension 128 × 128 pixels and grayscale values (8-bit color depth). The discrete model of the image degradation process under noise corruptions can be written as: whereū, f ∈ R N represent vectorized forms of the unknown clean image and of the observed corrupted image, respectively, while N ( · ) denotes the noise corruption operator, which in most cases is of random nature. In this work we considered two important types of noise, namely the additive (zero-mean) white Gaussian noise (AWGN) and the impulsive salt and pepper noise (SPN), which models saturated or dead pixels.
Denoting by Ω := {1, . . . , N } the set of all pixel positions in the images, for these two kinds of noise the general degradation model in (19) reads as AWGN : SPN : In case of SPN, only a subset Ω 1 of the pixels is corrupted by noise, whereas the complementary subset Ω 0 is noise-free. In particular, the corrupted pixels can take only the two possible extreme values V min /V max , where in our case assume 0 and 255 values, with the same probability. The amount of noise can be measured with an error rate computed as follows: ER % = number of corrupted pixels number of pixels in image × 100.
For what concerns AWGN, the additive corruptions n i ∈ G(σ, 0), i ∈ Ω, represent independent realizations from the same univariate Gaussian distribution with zero mean and standard deviation σ.
The performance has been evaluated by the Root-Mean-Square Error (RMSE) metric, defined as: whereū is the original reference image and u is the denoised output image. A lower RMSE indicates a more precise reconstruction.
For the remaining part of this work, when you come across the terms classical/quantum algorithm, we are referring to their respective implementation.
The selection of the regularization parameter λ, which exerts a crucial effect on the solution, has been carried out, for each test, by running the TV algorithm for a range of λ values in order to select by the trial-error strategy the optimal regularization parameter. Then the estimated selected optimal λ has been used in the quantum TV algorithm in order to compare the results of the two algorithms with the same optimal parameter.
A fundamental concept we must keep in mind when analyzing the obtained results, is that classical TV denoising uses floating point numbers, usually along with a value normalization, to be as accurate as possible. This leads to a more precise outcome and a finer quantization of the output image. On the other hand, our quantum computation uses integer numbers, so an approximation had to be applied (as previously described in Section 4). This difference has an impact on the output images, which present coarser improvements than the classical ones.
For testing purpose, we created an implementation of Quantum Median Filter using Qiskit, a Python library for quantum computing simulation [9]. The simulations of the quantum algorithm has ran on the Galileo100 supercomputer (CINECA), with the following cluster configuration:

Example 1: AWGN denoising
We consider the problem of denoising the three test images lena, QR, cameraman, corrupted only by AWG noise with standard deviation σ = {5, 10, 15}, as shown in the first column of Fig. 10, Fig.11, and Fig.12.
The denoised images obtained by applying TV and QTV algorithms are illustrated in Fig. 10, Fig.11, and Fig.12, column-wise for the three images. For each denoised image the RMSE obtained value is reported below.
From a visual inspection the denoised images from TV and QTV present a comparable quality, even if the RMSE highlights the lost in accuracy due to the mentioned integer arithmetic representation followed by QTV.

Noisy
TV QTV  Moreover, the quantum TV algorithm is able to denoise images corrupted by severe Gaussian noise, as illustrated in the last row of Fig. 10, Fig. 11, and Fig. 12 for σ = 15. Unlike, the quantum median filter proposed in [10], as any standard median filtering, is not appropriate to remove this kind of noise. The median filter is instead well-known to be an excellent image denoiser in case of salt-and-pepper noise because it does not blur the image, as a mean filter would do. However, despite its name, the median filter is not a filter because it does not respect the linearity property.

Example 2: SPN denoising
In this example we applied TV and QTV for the denoising of the three test images lena, QR, cameraman, corrupted by SPN noise with error rate ER % = {5, 10, 30}. The noisy images are shown in the first column of Fig.13, Fig. 14, and Fig.  15, together with the denoised images in the second (TV) and third (QTV) columns, along with the associated RMSE values, reported in the bottom.
From these results, we can see how well quantum algorithm performs when compared to its variational counterpart. QTV results present excellent qualitative performance, with minimal RMSE differences.
However, by a visual inspection of the denoised images, we notice that, for both the algorithms, some pixel clusters were not completely denoised. This can be due to either to the limited pixel neighborhood considered (4 pixels), or to the L 2 -norm metric fidelity in our model (5), when it is well known that SPN can be better treated by a L 1 -norm fidelity, which, anyway, leads to a non-differentiable fidelity term.

Conclusions and discussion
In this work a quantum approach is proposed for total variation denoising, and its corresponding quantum circuit is designed. Specifically the quantum TV implements the anisotropic median formula presented in [13]. The main idea of the approach is that first the classical image is converted into a quantum version based on the quantum representation (NEQR) of digital images, and then three quantum modules are applied to realize the neighboring collection for each pixel in the image, weight calculation, and median extraction. Finally, an image measurement process collapses the quantum state into a resulting denoised image. From the complexity analysis in Section 5 we derived a polynomial complexity in the number of pixels N for the NP module, and a polynomial complexity in logarithm of the number of colors N q = 2 q for the PC and MF modules. The experimental results show that the quantum TV performance is comparable to the classical variational TV approach. However, we highlighted several issues that need to be addressed to make the proposal a competitive QIP algorithm, as its variational counterpart. For example the Neighborhood collection module is an expensive operand due to the NEQR image representation. Even though NEQR is still one of the most used representation methods in QIP and the most suitable choice for this work, future developments should definitely search for other QIR alternatives, or develop more efficient versions of the same representation, such as parametric quantum circuits that take advantage of data structure for improving image processing activities.

Appendix A Quantum Gates for QTV
In this appendix we illustrate a few modules composed of basic quantum gates used in the design of QTV algorithm. We used classical logic for designing a quantum version of a Half-Adder (HA) and a Full-Adder (FA), as illustrated in Fig. A1(b) and A1(a) . These subcircuits need three additional auxiliary qubits, for storing temporary results needed for summing. FA does not provide a minimum/maximum cap over result: to solve this, we can use the temporary result stored in auxiliary qubits to manually fix the output to the correct value, as illustrated in Fig.A1(c).
The Comparator module is illustrated in Fig.A1(d). This module uses a very small amount of elementary gates and a reduced number of auxiliary qubits: if Reset option (RES) is not available, only q ancillas are needed for Comparator to work.