1 Introduction

Image segmentation [26] is one of the most important image processing tools, as it is frequently applied in an early stage of image analysis in order to extract some high-level features of the input image and significantly simplify further computations. Image segmentation is a process of dividing entire image into multiple separate segments – areas inside which pixels expose similar characteristics (e.g. intensity, hue, texture). It might be also defined as an activity of assigning each pixel of an image with a label, in a way that pixels sharing similar traits of interest are labeled alike. Image segmentation has found a broad range of applications: from a very trendy in recent years computer vision (object detection and recognition) [13] to medical diagnostics and treatment planning [14]. Despite many efforts and development of plethora of various methods, image segmentation remains an open research area, as there is still a need for inventing more efficient and more precise solutions.

A very interesting method for segmentation was invented by Grady [9]. The author proposes solution utilizing classical random walks, a very successful computational framework widely used for simulation of some natural and social processes in physics, biology and even economics [22]. The mechanism of random walk action is quite simple: a particle (walker) explores given space of possible states by moving at each time step to a new state according to current state and given transition probabilities. The walker behavior can be useful for determining some properties of the explored space. Based on the random walk there has been developed a analogous quantum model – quantum random walks. It has been shown that quantum walks can benefit from the extraordinary properties of nanoscale world (quantum superposition, interference or entanglement) to achieve quadratic or even exponential [7] speedup over their classical counterparts for some particular problems. Therefore it seems reasonable to explore this branch of quantum computation in the quest for new, better algorithms.

The aim of this research is to elaborate and examine a new algorithm that applies concept of quantum walks to image segmentation. It is inspired by results from its classical counterpart [9]. In this paper, we present concept, implementation and results of two versions of the solution. The results were obtained by simulation on classical computer, but the algorithms were designed in a way that will allow them to be run on a real quantum computer, when the actual device will be available. This approach seems to be promising, especially in a context of recent work that describes implementation of simple quantum walk on IBM-Q hardware [4].

The rest of the document is organized as follows: Sect. 2 discusses research concerning image segmentation and quantum walks. Section 3 presents elaborated algorithms, while Sect. 5 examines their performance and accuracy. Finally, the effort put in this paper is summarized in Sect. 6.

2 Related Works

Image segmentation has been a carefully studied problem over the years and there has been developed a multitude of different solutions [21, 26]. One of the simplest are threshold-based algorithms, which assign a label to each pixel according to its intensity value in comparison to one or more thresholds, e.g. frequently used Otsu method [19]. Another approach that has been successfully applied is based on clustering; it determines segments by aggregating pixels in groups (clusters) according to chosen metric [17]. Deep learning and neural networks, especially convolutional neural networks, also proved to be very successful in various image processing tasks, one of which is image segmentation [6].

The interesting class are graph partitioning methods. To this category belong methods like: normalized cuts [23], intelligent scissors [18] and isoperimetric partitioning [10]. The aforementioned Grady algorithm also belongs to the last category of image segmentation solutions. As the input, the user defines (apart from an image) a set of seeds – pixels which already have assigned labels corresponding to the image segment that they belong to. The algorithm then tries to determine for each unmarked pixel, what are the probabilities of reaching each seed as the first by a random walk staring from given pixel. Each pixel is then assigned with a label of the seed that is most probable to be visited first by the walker.

With the development of quantum computing there emerged a new field of research – quantum image processing. In [24] there were presented methods for storing and retrieving as well as segmenting images in a quantum system. However, the proposed algorithm requires specifying objects in image through as set of their vertices, which rarely is a case and obtaining those vertices is a complicated task itself. Two other methods for quantum image segmentations were described in [5, 12], but these are the threshold-based solutions which do not take into account the spatial locations of pixels, which results in rather poor outcome.

Since Ashwin and Ashvin [3] showed that quantum walks on line spread quadratically faster than their classical analogues, the quantum model has been studied to find algorithms that would allow for similar speedup. The most spectacular achievement in this field was presented by Childs [7], where the author achieves exponential speedup by applying a continuous time quantum walk in a graph formed of two binary trees. There have been presented several search algorithms that can provide quadratic speedup benefiting from quantum walks e.g. works by Ambainis et al. [2] and Magniez et al. [15].

Over the recent years, there have been developed many various algorithms utilizing quantum walks. Most of the works consider rather theoretical problems, therefore it might be beneficial to elaborate an algorithm that would directly solve a practical, commonly encountered task.

3 Walk Model

In this paper we propose two versions of the algorithm for image segmentation utilizing continuous time quantum walks. In continuous time quantum walks (CTQW) transition from the current state to the next one does not take place in discrete steps, but can occur at any moment, according to a \(\gamma \) parameter that describes the rate of spread of the walk.

As the input, the elaborated algorithms expect following data:

  • Image – two dimensional image of size \(M \times N\) with each pixel \(p_{ij}\) described by two indices: \(i \in \mathbb {Z}_M\) (index of the row) and \(j \in \mathbb {Z}_N\) (index of the column). Each pixel can be represented as a single value of luminescence intensity or a vector of intensities in separate channels (e.g. vector of length 3 for the RGB color model).

  • Labels \(L = \{1, 2, 3, ..., l\}\) – set of labels, one for each segment that the image should be divided into.

  • Seeds \(S = \{(p_1, l_1), (p_2, l_2), ..., (p_s, l_s)\}\) – set of pixels \(p_k\) with already assigned label \(l_k\) from the set L; there should be multiple seeds for each label.

  • Parameter \(\beta \) – a parameter responsible for highlighting the stronger differences between pixels; it is used while transforming the image into a graph.

  • Parameter \(\gamma \) – rate of spread of a continuous time quantum walk, used to construct the evolution operator.

  • Parameter T – duration of the quantum walk (number of applications of evolution operator).

Firstly, the input image is converted into a weighted graph with vertices corresponding to pixels of the image. From each vertex there are drawn four edges to the neighboring pixels. Finally, each edge is assigned with a weight describing the level of similarity between adjacent pixels \(p_{ij}\) and \(p_{kl}\). We use the formula proposed in [9]:

$$\begin{aligned} \begin{aligned} w_{ij,kl}={\left\{ \begin{array}{ll} e^{-\frac{\beta d(p_{ij}, p_{kl})}{\max d(p_{ij}, p_{kl})}}, &{} \text {if }p_{ij}\text { and }p_{kl}\text { are connected,} \\ 0, &{} \text {otherwise,} \end{array}\right. } \end{aligned} \end{aligned}$$
(1)

where \(d(p_{ij}, p_{kl})\) is a metric of pixel similarity defined in our paper as follows:

$$\begin{aligned} d(p_{ij}, p_{kl}) = \sum _{c} (p_{ij}[c] - p_{kl}[c])^2, \end{aligned}$$
(2)

where c are the consecutive channels of pixels.

The main idea to tackle the problem of image segmentation is shown in the Fig. 1 illustrating quantum walks for two example labels: green and red. First, a separate continuous time walk for each label is started and the walker explores the image according to its content. Then, based on which walker was the most willing to explore given region of the image, assign appropriate label to each pixel, which results in dividing the image into segments. The details of constructing the quantum walk described in this paper are presented below.

Fig. 1.
figure 1

Visualization of the concept outline. The algorithm input is shown in (a) where circles symbolize pixels with their intensity and seeds are marked as two pixels with green and red ring. In (b), image is transformed into weighted graph based on intensities of neighboring pixels (weights are denoted by the edge width). In (c), a walk is started from a superposition of seeds with given (here green) label (colorful dots denote the probability of measuring walker at given pixel) and, in (d), proceeds for several steps. In (e), when the walk ends the limiting distribution (CTQW-LD) or the last state (CTQW-OS) is calculated. In (f), similar procedure is repeated for each (here red) label. In (g), final distributions of each walk are compared and each pixel is assigned with label of the walk that had the highest probability of being measured at given pixel (h). (Color figure online)

Full size image

Position Space. The position space for the walk is a Hilbert space of size \(M \times N\); each basis state corresponds to a pixel of the image: .

Initial State. The initial state of the walker is a superposition of positions of all the seeds for given label. Let the current label be \(l_i \in L\) and the \(S_{l_i} = \{(p_1, l_i), (p_2, l_i), ..., (p_k, l_i)\} \subset S\) be a subset of seeds corresponding to that label. Then the initial state for walk corresponding to the i-th label has the following form:

(3)

Evolution Operator. The most important element of the quantum walk is its evolution operator, that determines the behavior of the walk. The current task is to construct an operator that would encourage walker to explore image according to the weights between neighboring pixel. The desired behavior is to make the walker advance preferably over the higher weights between similar pixels and avoid crossing big gaps in pixel intensities.

The general form of the evolution operator of continuous time quantum walk has the following form (here i means imaginary unit):

$$\begin{aligned} U(t)= e^{-iHt}, \end{aligned}$$
(4)

where t denotes time and H is a Hamiltonian matrix. The task is to find an appropriate Hamiltonian matrix H. It has to be a square matrix of size \(MN \times MN\) that would be hermitian (this ensures that the operator U(t) is unitary). It also has to yield satisfying results in terms of position space exploration by the walker according to the content of the image.

The idea is to construct the Hamiltonian based on the weights matrix defined by the Eq. 1 and a free real parameter \(\gamma \) that determines the rate of spread of the quantum walk:

$$\begin{aligned} H_{ij, kl}={\left\{ \begin{array}{ll} -\gamma w_{ij, kl}, &{} \text {if }i \ne k \vee j \ne l.\\ \gamma \displaystyle \sum _{k',l'} w_{ij, k'l'}, &{} \text {if }i = k \wedge j = l. \end{array}\right. } \end{aligned}$$
(5)

Please notice, that since weight matrix is symmetrical (\(w_{ij, kl} = w_{kl, ij}\)) and its values are real numbers, the constructed matrix H is a hermitian matrix and therefore the operator U(t) is unitary.

4 Algorithms Description

Based on the elaborated continuous time walk model, here are proposed two versions of the algorithm for image segmentation: with limiting distribution (CTQW-LD) and one shot version (CTQW-OS).

Limiting Distribution Version. Due to the application of unitary operators, quantum random walks do not converge to a stationary distribution, as it is the case with the classical random walks. Instead, it was proposed by Aharonov et al. [1] to consider average distribution and limiting distribution.

Definition 1

[1] Consider the probability distribution on the nodes of the graph after t steps of quantum walk starting from the initial state :

(6)

Then the average distribution is the mean over the distributions in each time step until T:

(7)

Notice that the average distribution can be understood as a measure of how long the quantum walker spends at each of the position states throughout the first T steps of the walk. It was shown in [1] that, for any initial state , the average distribution converges as the time approaches infinity and the limit is denoted as limiting distribution.

Definition 2

[1] Limiting distribution of the quantum random walk starting from the initial state :

(8)

The limiting distribution is not unique, but depends on the initial state.

figure a

The version of the algorithm using limiting distribution (CTQW-LD) is shown in the Algorithm 1. For each label the separate quantum walk is performed. After each step of the walk the probability distribution of measuring the walker at given position is recorded. In our simulations we could directly read the walker state after each step, as the calculations were performed using matrix algebra. Also after introducing the optimization (described below) we obtained the state after given number of steps by changing the equation parameters. However, on a quantum device retrieval of the quantum state is not straightforward, This is achieved by performing so called state tomography repeating the walk for given number of steps multiple times and making a measurement which results in finding the walker at some position on the image. The number of samplings grows with the number of basis states (the image size) and the desired accuracy of the reconstructed state. This allows to retrieve the probability distribution of the walker after each time step. Upon finishing each walk (the walk for each label) the limiting distribution of the walk is calculated (due to the fact that there is no possibility of running a walk for infinite number of steps, the limiting distribution is approximated with average distribution of T steps). Finally, the result is obtained classically, as each pixel is assigned with label of the seed, for which walker starting from that seed expressed the highest average probability of being measured at that pixel.

figure b

One shot version of the algorithm (CTQW-OS) is shown in the Algorithm 2. This solution is almost identical to the CTQW-LD, but instead of calculating the limiting distribution of each walk, only the distribution after last state is used. As it will be shown in the next section, this method gives results with a bit worse segmentation accuracy, but for careful choice of parameters the difference is ommitable. However, it does not require obtaining the probability distribution after each step, which is a tremendous profit, especially if the algorithm would be executed on a real quantum device. Quantum state tomography, which is used to retrieve the quantum state of the walker is an expensive operation, but in the CTQW-OS version of the algorithm it needs to be performed only once for each label, what is a great advantage over the CTQW-LD algorithm.

Optimization. In order to improve the performance of simulations of the algorithms on classical computers there have been introduced some optimizations. First, each walk was performed in parallel in separate processes. Additionally, it was tested that sampling the walker distribution not after each step, but for about every \(100^{th}\) step in the range [1, T] gives almost undisturbed results. Therefore, instead repeating evolution operator U, we used the numerical solver of Lindblad master equation [25].

5 Evaluation and Results

The evaluation of the proposed solutions could not, obviously, be performed on the real quantum device, as the requirements e.g. the number of qubits needed to express the walker state, the level of customization the quantum operators etc. are beyond the capabilities of currently accessible quantum computers e.g. IBM-Q [11] or Rigetti [20]. Therefore, there had to be performed a simulation of quantum computations on a classical machine. The project code has been prepared using Python language (version 3.4.3) with QuTiP open-source Python library for simulation of quantum systems (version 4.2.0). The program was executed on a computer with processor Intel Core i5-4210M and 16 GB of RAM. We used the Berkeley Segmentation Dataset and Benchmark (BSDS500) [16] dataset which contains a multitude of photographs, with corresponding contour images performed manually by people.

Fig. 2.
figure 2

Comparison of results on an example of color and grayscale test images. The original images are shown in (a, h) and the manual segmentations in (b, i). Results of Grady’s algorithm is shown in (c, j). Results of CTQW-LD with configuration: \(\beta =150.0, \gamma =0.001, T=5000\) is shown in (d, k), CTQW-LD with the most optimal configuration: \(\beta =150.0, \gamma =0.001, T=20000\) in (e, l), CTQW-OS with the most optimal configuration: \(\beta =150.0, \gamma =0.001, T=5000\) in (f, m) and finally, CTQW-OS with configuration: \(\beta =150.0, \gamma =0.001, T=20000\) in (g, n).

Full size image

The most important parameter for the proposed solution is the set of seeds that needs to be provided separately. We considered one set created manually, by careful choice of evenly spread pixels with a bit higher concentration around borders, especially the weak and unclear boundaries (these sets contained about \(0.2\%\) pixels of the whole image); and three sets that were prepared automatically by drawing uniformly a fraction of pixels from the image and assigning them with the labels based on the ground truth segmentation from BSDS500 dataset (with the concentrations at the level of \(0.2\%\), \(0.5\%\) and \(1\%\), respectively). There are also other three free parameters:

Parameter \(\beta \) – it is used during weights calculation; the greater its value, the more significant the stronger differences between pixels are. So the high \(\beta \) parameter makes the walker stay closer to its seed rather than disperse and explore further regions. However, the change of this parameter does not significantly influence the outcome segmentation.

Parameter \(\gamma \) – the rate of spread of the walker. The greater the value of this parameter, the greater is the area explored by the walker upon a single application of evolution operator. This means that one could set \(\gamma \) to a high value to limit the number of steps required. Indeed, increase of this parameter results in obtaining better outcome for smaller number of steps. But for relatively high values of \(\gamma \), as the number of steps grows, the segmentation accuracy drops almost as rapidly as it reaches its peak value.

Parameter T – determines the number of steps (number of applications of the evolution operator). Too small values result in poor image exploration by walkers and appearance of some blank areas. Too high values result in walkers diffusion to neighboring segments that reduces the segmentation quality over time. It is especially important for the CTQW-OS algorithm that strongly suffers from too high T value.

The optimal configurations for both of algorithms have been chosen as follows: \(\beta =150.0, \gamma =0.001, T=20000\) for CTQW-LD solution and \(\beta =150.0, \gamma =0.001, T=5000\) for CTQW-OS method.

In addition to the ground truth segmentation obtained from the manually created image outlines from BSDS500 dataset, we have also prepared a reference implementation of the classical random walks algorithm [9]. Each segmentation method was then compared against the ground truth images and finally there was calculated the percentage of pixels with consentaneous labels.

The representative results of proposed algorithms are presented in the Fig. 2, which shows segmentations of two sample images: a seashore landscape in color in the Fig. 2a and a man under a tree in grayscale in the Fig. 2h. Then there are reference segmentations: the Fig. 2b and the Fig. 2i are ground truth images included in BSDS500 dataset, while the Fig. 2c and the Fig. 2j were prepared using Grady’s algorithm. The last four images for each photo were obtained with algorithms proposed in this paper: the Fig. 2d and k present results of CTQW-LD solution with configuration: \(\beta =150.0, \gamma =0.001, T=5000\), Fig. 2e and l come from CTQW-LD algorithm with optimal configuration: \(\beta =150.0, \gamma =0.001, T=20000\), while Fig. 2f. Figure 2m were obtained with CTQW-OS method with optimal configuration: \(\beta =150.0, \gamma =0.001, T=5000\). Figure 2g and n show CTQW-OS with configuration: \(\beta =150.0, \gamma =0.001, T=20000\). As it can be seen, 5000 steps is too little for the CTQW-LD algorithm (visible non-covered black areas in Fig. 2e and l). For CTQW-OS, if the walk duration is too long, the walks from adjacent segments tend to diffuse causing deterioration of the segmentation quality (Fig. 2g, n). As could be expected CTQW-OS is a bit inferior to CTQW-LD. Nevertheless, the difference is not obvious, due to careful choice of parameters (especially number of steps), while the gain in terms of performance is tremendous. CTQW-LD looks comparably, if not better than Grady’s method (please notice the difference around the head of the man under tree). Obviously, all of these solutions are imperfect which can be clearly seen in the areas of weak and blurred margins (like, the line of horizon or leaves of the tree).

Table 1. Average accuracy (over test images) of developed algorithms as well as Grady’s method for different sets of seeds.
Full size table

Table 1 presents average accuracy (over test images) of developed algorithms in comparison to Grady’s method for different sets of seeds. It clearly shows that the segmentation accuracy increases with the size of seeds set. Manually chosen seeds give better results than the same amount of automatically drawn ones, due to their careful selection. More accurate of the proposed solutions is CTQW-LD algorithm, which gives almost as good results as reference algorithm or even outperforms it in some cases. The other method provides slightly worse accuracy. The last two rows show how sensitive to the walk duration is the CTQW-OS solution; the accuracy drops significantly if the number of steps is too high. Nonetheless, for optimal configuration of parameters the results of both developed methods might be considered satisfactory. The only unfavorable aspect is the quite considerable amount of seeds needed to obtain a decent segmentation accuracy. As the \(0.2\%\) seems to be a bit too little, there has to be specified about \(0.5\%\) seeds, which for image of size 320 \(\times \) 480 means that over seven hundred of pixels need to be labeled.

The average execution times for test images of size 320 \(\times \) 480 and optimal configurations were as follows: about 5–6 min for CTQW-LD and around a minute in case of CTQW-OS method. Along with the growth of \(\gamma \) or T the duration of the calculations increases linearly. For comparison Grady algorithm for the same samples requires about 30 s. This means that developed solutions are comparable both in terms of accuracy as well as performance to Grady algorithm. Of course, these are the execution times of simulation on classical computer which would not practically matter during transition to a quantum device, nevertheless the proposed results, especially CTQW-OS, are expected to do well also in quantum conditions.

6 Summary, Conclusion and Future Work

In this paper, there have been presented two novel quantum algorithms for performing image segmentation based on continuous time quantum walks: one calculating the limiting distribution of the walk and the other utilizing only its last state. Both methods have been closely examined and tested using simulation of quantum computations on a classical computer. The obtained results are satisfactory, of the same accuracy and similar performance as the reference classical Grady’s [9] algorithm. A crucial difference between simulation and the actual execution of the algorithm on a quantum computer is the fact, that during simulation, there is always access to the whole quantum state. In real conditions the actual state of the particle is not known until measurement, which reduces it to one of the basis states according to probability distribution determined by the current state. So, to gain the knowledge on the walker state after a given number of steps, one needs to perform the same walk multiple times and measure it in order to reconstruct the probability distribution. In these conditions the CTQW-LD solution is not very effective in quantum case – in order to obtain the limiting distribution the experiment has to be repeated multiple times for different numbers of steps. For this reason, it would be vastly beneficial to apply the CTQW-OS algorithm, which requires the construction of the probability distribution only once per walk.

The segmentation presented in this work is a seeded-based method. User has to specify a set of pixels with assigned segment to which they belong which is not an automatic approach. Also to achieve good results a considerable amount of seeds needs to be provided. Therefore, future work in this field might concern elaboration of either a way to limit the number of required seeds or a method for automatic specifying the seeds, similar to those currently applied in medical diagnostic images segmentation [8]. This would allow to significantly increase the user experience while using the proposed solutions for image segmentation.

The future work on this matter could consider translating the algorithm (for some small size images) on a quantum device - firstly on a quantum simulator like IBM Q simulator or Atos QLM and further to make an attempt to run it on a quantum computer (e.g. IBM Q or Rigetti). This would require disassembling the evolution operator into basic one and two qubit quantum gates, as only those are handled by current quantum devices, as well as taking into account error correction strategies to mitigate the risk of quantum state decoherence.