Synthetic Aperture Radar Image Segmentation with Quantum Annealing

In image processing, image segmentation is the process of partitioning a digital image into multiple image segment. Among state-of-the-art methods, Markov Random Fields (MRF) can be used to model dependencies between pixels, and achieve a segmentation by minimizing an associated cost function. Currently, finding the optimal set of segments for a given image modeled as a MRF appears to be NP-hard. In this paper, we aim to take advantage of the exponential scalability of quantum computing to speed up the segmentation of Synthetic Aperture Radar images. For that purpose, we propose an hybrid quantum annealing classical optimization Expectation Maximization algorithm to obtain optimal sets of segments. After proposing suitable formulations, we discuss the performances and the scalability of our approach on the D-Wave quantum computer. We also propose a short study of optimal computation parameters to enlighten the limits and potential of the adiabatic quantum computation to solve large instances of combinatorial optimization problems.


Introduction
Synthetic Aperture Radar (SAR) is an imaging technique that utilizes radar signals to create highresolution, three-dimensional reconstructions of terrain and objects by exploiting the motion of the radar antenna and signal processing algorithms [1].Images obtained with SARs cover dozens of kilometers with a resolution between 15 and 30cm, resulting in images of billions of pixels.Due to the low interpretability and noisiness of SAR images, radar operators are using Automatic Target Detection (ATD) algorithms.The goal of these algorithms is to associate a group of pixels to an object of interest, as a vehicle or a radar station, by segmenting the image in different regions.In real time area surveillance, due to the size and the high number of images coming from multiple sensors, current computing capabilities and the restrained number of qualified radar operators greatly limits the number images that can be processed.Image segmentation is the process to divide an image in multiple classes, called segments or regions, in order to extract information about the image.Segmentation results takes the form of a set of labelled pixels, where each label corresponds to a class.The situation where the true label of each pixel is not known is sometimes called unsupervised image segmentation.In the literature, unsupervised image segmentation can be done with region smoothing [2], region growing [3], clustering [4,5], watershed segmentation [6], graph-based methods [7].Due to the exponentially increasing number of possible labellings, solving large instances of image segmentation problems often result in sub-optimal sets of labels.Quantum computing is an emerging technology which exploits the laws of quantum mechanics in order to perform logical operations.Instead of classical bits, quantum computers operate on qubits, which are in a superposition of two states.Among various applications, solving optimisation problems with quantum devices promises an exponential speedup over classical approaches [8], especially for large instances of NP-Hard problems.There are currently two main approaches in the design of quantum computers: Circuit-oriented quantum computers and quantum annealers.Circuit oriented quantum computers have a sequential approach of quantum computation, using gates to perform operations on single or multiple qubits.Quantum annealers have a simultaneous approach of quantum computation, making all the qubits involved in the computation converge from an initial state to a final state.In this work, we aim to propose an approach to speed up the processing of a high number of high resolution SAR images in order to tackle the current limitation regarding the analysing capability.Beyond this latter problem, our approach could also contribute for the training of Automatic Target Recognition (ATR) classifiers [9], by providing large datasets of segmented SAR images obtained in a reasonable time.Although having, today, more qubits than their circuit oriented analogs, quantum annealers are limited to optimisation problems in the form of Quadratic Unconstrained Binary Optimisation (QUBO) problem.In image segmentation, graph based methods consist in representing the image as a graph, where the nodes represent the pixels and the edges represent the neighbour between pixels.Hence, segmenting the image is equivalent to a graph problem, like coloring problems [10] or graph cut [11].Among these graph methods, Markov Random Fields (MRF) provide an effective way to model the spatial dependencies in image pixels [12].The goal of this paper is to pursue this work by proposing a QUBO formulation for MRF in order to perform unsupervised image segmentation, and compare the performances of our quantum approach with classical results and methods.We also propose an Expectation-Maximisation (EM) [13] inspired approach to achieve a satisfying unsupervised segmentation of Synthetic Aperture Radar (SAR) images.In the end of the paper, we propose a comparison of our algorithm with non-quantum approaches, and we discuss about its performance and scalability.

Markov Random Fields
MRF are probabilistic graphical models representing a joint probability distribution over a set of random variables.In an MRF, each node represents a random variable, and the edges between nodes represent statistical conditional dependencies between variables.Let X be a set of observables and Z a set of hidden variables on which one aims to infer.Following the Bayes' inference rule [14], we have : with θ a given set of parameters characterizing the distribution.In the following, we denote by x = {x i } i∈{1,...,n} the gray scale intensity of a pixel and z = {z i } i∈{1,...,n} its label.The associated random variables counterparts are denoted by capital letters X and Z respectively.The Maximum A-Posteriori (MAP) p(Z = z | X = x) (noted p(z | x) in the following for the sake of simplicity) is derived as follows: By identification from (2), we pose ψ i (z i ; θ) and ϕ i,j (z i , z j ).For ϕ i,j (z i , z j ), we make the hypothesis that the joint probability distribution p(z) can be decomposed in a product over the pairs of neighbouring nodes [15], which leads to the following general formulation: We note V i the ensemble of neighbouring nodes of the node i.The left term concerns the compatibility of the observable x with the hidden variable z.It is the "unary" term of the MRF corresponding to the log-likelihood of node i to be associated with hidden variable z i knowing its observable x i .The second term is the "pairwise" term of the MRF, which characterizes the compatibility of neighbouring hidden variables.Hence, making the MAP estimation of the MRF is equivalent to find the set of hidden variables z minimizing the sum of both terms, which is NP-Hard [16].In the literature, classical algorithms as belief propagation [17] or variational methods [18] are used to find the optimal sequences of hidden variables.In practice, for a large number of observables, computation time explodes as the number of possible hidden variables sequences exponentially increases [16].

MRF for Image Segmentation
In the following section, we will consider adapting the MRF model presented above for segmenting a grey-scale image of N pixels in Q distinct classes.The observable x i corresponds to the intensity of the i th pixel and its label is z i .Pixel intensity takes values from 0 (black) to 255 (white).Our input data is the set of pixels i ∈ {1, . . ., N }, the set of labels q ∈ {1, . . ., Q} and the set of neighbouring pixels V i for a pixel i. Let's pose the energy function H(z) such as : The goal is to find a parameterization of the functions ψ i (z i ) and ϕ i,j (z i , z j ) to make the minimum of H(z) correspond to a satisfying segmentation of the image.Here, ψ i (z i ) can be interpreted as the function measuring the cost of assigning the label z i to pixel i knowing its intensity x i .ϕ i,j (z i , z j ) can be interpreted as a function measuring the cost of assigning labels z i and z j to the two neighbouring pixels i and j.From (2) we deduce : with q ∈ {1, . . ., Q} and I i ∈ {0, . . ., 255} the intensity of the pixel i.The minimal value of this term is attained for p(x i | z i ) maximal.Hence, this term favors that each pixel is associated to the class he most likely belongs to by only considering its intensity.From (4), we pose the following pairwise term: with B i,j a real positive number and δ(z i , z j ) = 0 if and only if z i = z j , δ(z i , z j ) = 1 else.Literature provides different values for B i,j .In the Potts model [19], B i,j is defined as a constant value.In other models like Cauchy [20] or Huber [21] models, B i,j is a function of the intensity of the intensities of x i and x j .In this paper, for the sake of simplicity, we will consider the Potts model.In the result section, we will discuss the setting of B i,j This intuition behind this formulation comes from the hypothesis that connected pixels are more likely to have similar labels in a neighbourhood.Hence, the pairwise term favors the fact that all the pixels of the image have the same label, but combined with the unary term, it prevents the labelling of isolated pixels in the wrong class.

Expectation Maximisation Algorithm
The Expectation Maximisation (EM) algorithm is an unsupervised iterative algorithm used to make a Maximum Likelihood Estimation (MLE) [13] of statistical models parameters with missing data.In our work, we make the hypothesis that, for each region, the distribution of x i follows some known probability distribution model, and use EM to estimate the probability distributions for the unary term of the MRF.The EM algorithm can be decomposed in two steps : the expectation step (E-step) and the maximisation step (M-step).In the E-step, the algorithm computes the cost function of the problem for a given set of parameters θ t .This step returns a set a hidden variables and its corresponding cost i.e. energy.The M-step consists in updating the problem parameters, by maximizing the expected log-likelihood computed during the E-step.By alternatively repeating these two steps, the algorithm converges to a local maximum of the likelihood function.Even if demonstrations of the theoretical convergence of the EM algorithm have been achieved [22], EM remains sensitive to the value of θ 0 i.e. the value of θ t at step 0, which may lead to sub-optimal solutions [23].However, in our problem, as we do not have the exact ground truth for SAR images, our requirements of quality for the solutions allow us to be satisfied by near-optimal values of the estimated parameters θ.

Problem Implementation
As we have seen in section 2.2, we have on one hand MRF which are powerful models for image processing, embedding both unary and pairwise interactions between pixels.On the other hand, EM is a algorithm used for parameter estimation of statistical models, but with a computationally expensive Expectation step (E-step).Previous work shows combination of EM and MRF to obtain satisfying segmentation.[24] proposes an approach where problem parameters are estimated without consideration of their spatial distribution, then executes a MRF to infer the pairwise interaction term.[25] includes parameters based on the texture of images in addition to the color/intensity observables.In the following sections, we propose two QUBO formulations for the MRF image segmentation.The first one is a two-classes segmentation approach and the second one is a generalized version for Q classes.

QUBO formulation for 2-classes segmentation
Is this section, we propose a QUBO formulation to segment an image in two regions, formulated as a graph cut problem.Graph cut problems are a class of optimization problem aiming to partition a graph into different subsets, such that some cost is minimized [11].In the literature, algorithms as Grab Cut [26] or Alpha-Expansion [27] are used to segment images through graph cut methods without any quantum hardware implementation.In this paper, a cut is characterised by two neighbouring pixels i and j having two different labels i.e. with z i ̸ = z j .If the cut results in a increase of the cost function value, we speak of penalty.If the cut results in a decrease of the cost function value, we speak of a bonus.Here, we consider an image of N pixels.For a pixel i, we note V i the set of its neighbouring pixels.To each pixel is associated a label z i ∈ {0, 1} and a grey scale intensity x i ∈ {0, 255}.Label value z i = 1 implies that the pixel belongs to the "object" class, and z i = 0 implies that it belongs to the "background" class.For each value of x i , we associate the posterior probability p(x i | z i = 1) and p(x i | z i = 0) such as p(x i | z i = 0) + p(x i | z i = 1) = 1.If two pixels are neighbours, there exists an edge between their two respective nodes in the graph.
In order to implement the unary term, we also introduce two ancillary pixels a and b with respective labels z a = 1 and z b = 0. We consider that all pixels of the image are neighbours to these two ancillary pixels.Thereafter, in the graphical expression, we will call the nodes associated to these pixels the "ancillary" nodes.There exists an edge between each ancillary node and image nodes.
Following this formulation, each image node is bound to the nodes of its neighbours in the image plus to the two ancillary nodes.Our goal here is to parameterize the cost of cutting these edges, resulting in the global cost function of the graph cut problem, in the form of QUBOs.In the following sections, we will refer to the different cost function as Hamiltonians, in reference to the quantity defining the energy of a system in quantum mechanics First, we need to guarantee that z a = 1 and z b = 0. Hence, we pose the linear "constraint" Hamiltonian: h A (z a , z b ) is minimal if z a = 1 and z b = 0. We also define the "cut" Hamiltonian as follows : This Hamiltonian is equal to 0 if z i = z j , and equal to 1 if z i ̸ = z j .Then, we define the "unary" Hamiltonian as follows : If Hence, the pixel is labelled as the class it has the highest probability to belong.The global minimum of this Hamiltonian corresponds to each pixel being labelled to the class it has the highest probability to belong, which satisfies the formulation of 2.2.For the pairwise term, we pose the following "pairwise" Hamiltonian : In this expression, a penalty of value 1 is applied if z i ̸ = z j .From ( 7), ( 9) and ( 10), we pose the "problem" Hamiltonian as follows : With λ P and λ A two positive real multipliers.From (6), we pose λ P = B a positive real constant.For λ A , we have to ensure that it is always favorable, in term of value of the cost function, to respect that z a = 1 and z b = 0 for all labelling z.To do so, let's consider the extreme scenario where all the pixels of the image are labelled in the "object" (resp."background") class, and all pixel have an equiprobability of belonging to class "object" or "background".Let's suppose z a = 1 and z b = 0. Following the above conditions, ∀i, z i = z a and p θ (x i | z a ) = p θ (x i | z b ) = 0.5, Hence, we have here the maximum possible minimal value of (9).In that case, setting z b = 1 would apply a bonus of N log(0.5) to the cost function.In order to guarantee that this configuration will never be favored, we have to set λ A > −N log(0.5).Hence, ∀z, with z a = z b , H U (z, z a , z b ; θ) + λ A h A (z a , z b ) > 0, which ensures that any sequence minimizing H(z, z a , z b ; θ) also minimizes h A (z a , z b ).

QUBO formulation for Q-classes segmentation
In this part, we extend the formulation presented in 3.2 for Q classes.Lets consider an image of N pixels that we want to divide into Q disjoint subsets.In order to encode the Q possible labellings for each pixel, we pose the vector z ′ ∈ {0, 1} N Q , and we note z ′ i each bloc of length Q that is assumed to have only one non-null value.In the following, we define the function ψ(i, q) = (i − 1)Q + (q − 1) such as z ψ(i,q) is the i th bloc of z having its non-null value at the q th index.This method is called one-hot encoding, and we will consider that z ψ(i,q) = 1 implies that pixel i is labelled q.The rest of notations remains unchanged.
In order to ensure one-hot encoding for all, we define the "one-hot" Hamiltonian as follows : As in the binary case, we have to define a single ancillary pixel α such as z α = 1.To ensure its value, we pose the following "constraint" Hamiltonian : The "unary" Hamiltonian for the Q-classes is defined as follows : Following a similar principle as the unary Hamiltonian of 3.2, the goal is to maximize the sum of bonuses, by cutting the most negative edges.Here, one-hot encoding imposes that there are Q − 1 cuts between nodes corresponding z ψ(i,q) ∀q and z α for a given i.The optimal cut is achieved when the Q − 1 most negative edges are cut, which corresponds to the Q − 1 lowest p θ (x i | z ψ(i,q) ).Hence, the only non-cut occur for the maximal value of p θ (x i | z ψ(i,q) ) which corresponds to maximum likelihood term in 2.2.
The "pairwise" Hamiltonian for the Q-classes is defined as follows : Considering that one-hot encoding is respected, if pixels i and j have the same class, then z ψ(i,q) = z ψ(j,q) ∀q and no penalty is applied.From ( 12), ( 13), ( 14) and ( 15), we deduce the global QUBO formulation for the Q-class segmentation: With λ P ′ , λ OH and λ A ′ three positive real Lagrange multipliers.For λ P ′ , as the pairwise interactions follow the same hypothesis, we set λ P ′ = λ P = B.For λ A ′ , by following the same reasoning as in 3.1, we can again consider an extreme scenario.In the Q-class state, an equiprobability of belonging to each class implies that p(z ψ(i,q) | x i ) = 1/Q ∀q.Because there are Q − 1 cuts per pixel, the total number of cuts between the image nodes and the ancillary node is N (Q − 1).Note that this number of cuts is independent from the labelling.From the results of 3.2, we deduce that we have to set λ A ′ > −N (Q − 1)log(1/Q) in order to guarantee that the constraint is always respected for an optimal value of H ′ .
For the minimal value of λ OH , we have to consider the situation where neighbouring pixels i and j have different labels p and q.In this situation, a penalty of B is applied, which corresponds to the sum of two cuts.The first cut is between z ψ(i,p) = 1 and z ψ(j,p) = 0 and the second between z ψ(i,q) = 0 and z ψ(j,q) = 1.In this configuration, in order to guarantee that it is never favorable to set z ψ(i,p) = z ψ(j,p) and/or z ψ(i,q) = z ψ(j,q) , we have to set λ OH > λ P ′ 2 .

EM inspired approach for parameter estimation
In this part, we address the parameter estimation problem of the unary term of the MRF.For the sake of generality, the formulation we propose is based on the Q-classes formulation of 3.2, as setting Q = 2 is equivalent to a binary class formulation of 3.1.In this section, we propose an Expectation Maximization algorithm [13] to estimate the parameters of each distribution, then generalize it to any distribution and infer a general formulation for (15).In this section, we describe the algorithm used to obtain an unsupervised segmentation of an image of N pixel in Q segment.Pixel intensities with one segment are assumed to be a sample of a normal distribution with parameters θ q = (σ q , µ q ).For each step t of the EM algorithm, we note θ t q the distribution parameters of segment q at step t and we pose θ t = {θ t 1 , . . ., θ t Q }.The initial distribution parameter of the algorithm is noted θ 0 .The literature underlines the importance of well initializing θ 0 q in order to efficiently converge to a global minimum [23].In our approach, we used a k-means algorithm [28] in order to initialize the values θ 0 q for Q clusters.This algorithm provides a value of µ 0 q for each class, and a single value of variance σ for all clusters.For the sake of simplicity and computation cost, we consider that setting σ 0 q = σ ∀q is a satisfying initial variance value for all segments.For the E-step, we compute an expectation of the MRF for the set of parameters θ t .As a result, we obtain Q subsets of pixels intensity values x t q , one for each segment.For the M-step, the goal is to find the set of parameters θ t+1 such as ( 16) attains new minimum.If p θ (x i | z i ) follows a Gaussian model, this minimization is explicit [29], otherwise, minimization has to be performed numerically [30].In order to evaluate if convergence is reached, we compute with L(θ t ) = H ′ (z, z a ; θ t ) i.e. the equation ( 16) parameterized by θ.If ∆ t is least than some predefined threshold δ > 0, we consider that we have reached convergence, an retains the set of labels z t obtained during the E-step.Else, t ← t + 1 and the algorithm loops back to a new E-step.

Adaptation for SAR image processing
In this part, we consider adapting the above formulations to SAR images.In the literature, previous work assimilates the distribution of values of pixels in SAR images segments to a Weibull distribution [31,32], which is a special case of the general gamma distribution (GGD).The GGD has the following probability density function : With k, p, λ and t positive real numbers.In the above formulation, k and p are the shape parameters, λ is the scale parameter and t is the location parameter.In the case of a 2 parameters Weibull distribution, k = p and t = 0. Straightforward calculations lead to the following distribution: Unlike the Gaussian distribution, there is no explicit method to estimate the parameters of a Weibull distribution.Hence, for the M-step of the MRF algorithm described in section 4.1, we will rely on basic numerical optimization to update our problems parameters for each class.

Experimental context
In this section, we propose and discuss the results obtained on the D-Wave adiabatic quantum computer [33].D-Wave is a Canadian company which proposes access to their quantum devices through a cloud based service called D-Wave Leap.The device we use in the following is the D-Wave Advantage quantum computer, having 5760 available qubits and 15 couplers per qubits.The number of couplers quantifies the number of interconnections available between physical qubits.In the QUBO formulations presented above, a product between two variables implies there is a coupler between their two corresponding qubits.As this number of products can exceed 15 in our formulations, the embedding function provided by D-Wave libraries binds qubits with each other, creating logical qubits from multiple physical qubits.In the computation, logical qubits act as a single variable.Consequently, due to this embedding constraint, the actual limitation of problem size that can be implemented on the quantum machine is below the theoretical limitation (N 2 + 2 < 5760 for example for the bi-classes formulation).Recent publications from D-Wave present hardware having 7800+ qubits available and more interconnections.However, in the following, we will consider current limitations of D-Wave hardware at 5760 qubits and 15 interconnections.In order to compare our results, we also used a simulated annealing library provided by D-Wave called d-wave neal.We use this tool in order to implement larger instances of the problem which cannot be implemented in the quantum hardware.More precisely, we chose to simulate our results through the d-wave neal routines in order to segment images of a sufficient size to highlight the performances of our approach.For 2 classes segmentation, the maximum size of instances we can implement on the quantum computer is (22,22) with the binary formulation and (15,15) for the multi class formulation (with Q = 2).For Q = 3, the maximum image size is (9,9).Through our experiments, it appeared that the major limitation in term of embedding was the maximal number of connections between variables.As we had no access to D-Wave's embedding routines in order to find better formulations to our problem regarding the hardware constraints, we did not focus on finding an optimal embedding in our work.In comparison, d-wave neal enabled us to obtain results for larger images, up to (100, 100) for 2 classes segmentation and (80, 80) for 3 classes segmentation, with a limitation in term of memory on our classical devices.
Unfortunately, as we write these lines, d-wave neal does not provide a quantum noise simulator, resulting in over-optimistic results.Consequently, our estimations of algorithm performances will be based on an extrapolation of performances obtained for small images.In order to put the results in perspective, the following sections propose a quantitative and qualitative comparison with a region smoothing algorithm presented in [2].In this paper, the authors compare their approach with other unsupervised segmentation algorithms for both gaussian noisy and SAR images.Their results appear to be significantly better in computation time and accuracy on toy examples compared to fuzzy approaches.We also tested the fuzzy clustering algorithm proposed in [5], which gave satisfying results for gaussian noisy images but did not show good results on SAR images.
Concerning the weight B of the pairwise term, search heuristics presented in Fig. 1 have shown that best values were obtained for λ P around 0.5, which corroborates the results presented in [19,34].Moreover, it appeared in our tests that fine tuning the value of λ P as a function of Q would improve the quality of segmentation.Hence, in all the following test, we set λ p = 0.6, λ p = 0.5 and λ p = 0.35 for resp.Q = 3, Q = 4 and Q = 5.In the following, segmentation accuracy is defined as the percentage of well-labelled pixels.The computation time for a single sample (the execution time on the quantum annealer) is the sum of the delay time, the annealing time and the readout time.Delay time is fixed at 21µs, and consist in initializing the system in an equiprobability of obtention of any sequence.Then, the annealing time represents the convergence time from the initial state to the final state.The latter is defined by the QUBO, and its minimal energy configurations correspond to solutions minimizing the criterion.Previous work shown that the annealing time logarithmically scales with the number of variables [35].In our work, for the sake of simplicity, we will consider an annealing time of 200µs regardless the number of variables.Concerning the readout time, it logarithmicaly scales with the number of variables, for a value in the order of 100µs.Due to the inherent error of quantum computation, we consider repeating the annealing, and then keep the best result i.e. the set of labels with the minimal energy among all trials.Hence, repeating the number of trials increases the chances to obtaining an optimal solution.In section 6.4, we will discuss the settings to increase the probability of obtention of optimal solutions whilst minimizing the hardware related error.

Results for Gaussian distributions
In this section, we consider segmenting images with an additional Gaussian noise.For binary and multi-classes segmentation, we will respectively use formulations ( 11) and ( 16) as the E-step of the EM algorithm.Due to the variable queue time to access the D-Wave hardware, we consider that the computing time of the E-step is equivalent to the computation time for a single sample times the number of samples.Hence, the total computation time is equivalent to the initialization time (problem embedding + k-means algorithm), plus the number of steps of the EM-algorithm times the computation time of the E-step and the M-step.In order to get a better result and limit the impact of error at the end of the algorithm, once convergence is reached, we re-compute the E-step for a greater number of samples, and keep the labelling corresponding to the minimum energy as our result.
Fig. 2 presents qualitative results for generated images with Gaussian noise using our approach and the region smoothing algorithm proposed in [2].In Table 1, quantitative results are shown for different generated checkerboards with Q = 3 to Q = 5 classes.It presents the improvement of accuracy provided by the EM algorithm, as the parameters of each Gaussian corresponds to the ones of the generated image.For Table 1, Table 2 and Fig. 2, we present results obtained with the d-wave neal sampler, as hardware embedding constraints (especially one-hot encoding) are not matched for a high number of pixels and sequences respecting the one-hot encoding constraints for all pixels are rarely obtained for large instances (nb. of logical variables > 250).Nevertheless, for the binary formulation, as no one-hot encoding is required, embedding errors results in a loss of accuracy, without violating the constraints.For all the test on d-wave neal, the number of samples to 100 samples, the convergence criteria to δ = 5, the annealing time to 250µs and the maximal number of epochs to T = 30.
Quantitative  1: Average accuracy for multi-class segmentation for images with added Gaussian noise.In the above, Q corresponds to the number of class and σ corresponds to the variance of normal distribution of the Gaussian noise.Nb. var.corresponds to the number of (logical) qubits required to implement the problem.T represents the number of iterations.
Quantitative results presented in Table 1 underline the advantage of using the EM Algorithm to improve the quality for the segmentation.It appears that this improvement is proportional to the number of classes, the size of the image and the value of σ.Regardless the embedding constraints, computing the E-step for 100 samples on the quantum annealer would take at most ≈ 100 × 200µs = 0.02s plus ≈ 0.2s to update the cost function and re-implement it on the quantum hardware.For the Mstep, the approximate computation time is 0.05s.These results leads to an approximate computation time of ≈ 0.3s per epoch, which can be even more reduced by a faster classical computation.This duration logarithmicaly scales with the number of variables, as the optimal annealing time and readout time logarithmically scales with the number of qubits.As a comparison, d-wave neal sampler has a computation time of 5.3 seconds per epoch (image of size (40,40) and Q = 4) for the E-step (cost function update + computation) for 100 samples, which is consistent with the results of [36,37] obtained with non-quantum approaches.
In Table 2, we compare our results with region smoothing and fuzzy clustering algorithms resp.proposed in [2] and [5].Simulated annealing results obtained with d-wave neal underline the computation cost of sampling a MRF with classical approach.The corresponding computation time for each instance would be ≈ 9 seconds (30 × 0.3 seconds, average computation time mentioned above), time that could be consequently reduced by a faster computation of non-quantum calculus.For each epoch, sampling the MRF only takes 0.02s of quantum computation for 100 samples, instead of more than 100s in simulated annealing.This implies that processing larger images would not significantly increase the computation time.Moreover, assuming that quantum annealing would provide the same results   [2] and [5].Simulated annealing results are obtained by running our algorithm on 30 epochs using d-wave neal for the E-step.Notation t corresponds to the computation time in seconds.
as simulated annealing, for the largest instances, the algorithm provides better results in terms of accuracy than the two other approaches.These latter assessments lead us to think that our algorithm would be relevant for segmenting large SAR images, as described in the introduction.

Results for SAR images
In this section, we present results on SAR images with the adapted parameterization (18) described in section 5.For these tests, we will only provide qualitative results, as ground truth is not available.First results are proposed on cropped images of the MSTAR dataset [38], composed of various military vehicles SAR images similar to Fig. 3.As evoked in the introduction, datasets composed of MSTARlike segmented images can be used to train ATRt classifiers [9].Hence, our goal is to subdivide these images in 3 segments, corresponding to the vehicle, its shadow and the rest of the image.The object corresponds to a connected area composed of high intensity pixels and the shadow corresponds to a connected low intensity pixels area next to the object.For the chosen MSTAR images, we consider that a satisfying segmentation is obtained if the object and shadow segments are both connected areas next to the center of the image, without including background labelled pixels.As we aim to segment the image into 3 classes, the E-step of the EM algorithm will be executed with the (16) formulation, adapted to SAR images.Concerning the initialization of the EM algorithm, k-means and Weibull mixture models did not provide an accurate set of parameters, which led to nonsatisfying segmentation.Hence, we chose to initialize our model with a thresholding method.Thus, we consider that pixels with intensity of x i ≤ 7 most likely belong to the "shadow" class and pixels with an intensity x i ≥ 20 most likely belong to the "target" class.For other intensity values, we consider that they most likely belong to the "background" class.
In Fig. 4, we present results obtained with the d-wave neal sampler on some MSTAR dataset images, and compare its performances with the region smoothing [2] and fuzzy clustering [5].In this figure, rightmost images, obtained with the fuzzy clustering algorithm, do not present a satisfying segmentation of SAR images, as the shadow region is associated to the background.On images obtained with the region smoothing algorithm, target and shadow regions are well located but are not connected and roundish, which may complicate the recognition of the target.In Fig. 5, we present the results obtained with the d-wave neal sampler on the SAR Ship Detection Dataset (SSDD) [39].As in the previous figure, we compare the performances of our approach with the region smoothing [2] and fuzzy clustering [5] algorithms.Segmentations b/c/d are obtained on a (100, 100) image (a) with two classes (corresponding to the ship and the sea).For the two other segmentations, we consider (75, 75) images (e and i).Those images are segmented in 3 classes ranked from the brightest to the darkest: the ship, the land and the sea.On these figures, we highlight the efficiency of our algorithm for removing noise and artifacts on lower-quality SAR images (in comparison to MSTAR dataset images) whilst distinguishing each class with a satisfying accuracy.As in Fig. 4, evaluated classical approach present less satisfying results, by being less accurate in the preservation of the shape of the ship and the delimitation between the land and the two other classes.
For the same reason as in section 6.2, embedding constraints and hardware error greatly limits the size of the instances we can implement on the quantum annealer.Due to the non relevance of segmenting small images (of maximum size (15,15)), we chose to only present our results obtained on d-wave neal as a realistic expectation of future quantum hardware performances, especially in terms of embedding capability and error tolerance.

Discussion on optimal parameters and limits of the model
In the previous sections, our results underlined the current limitation of quantum technologies regarding the error rate and the embedding constraints.Nevertheless, we can already discuss on optimal parameters for the model as the value of the weighting of the constraints, and also optimal parameter for the computation, as the annealing time.In Fig. 6, we can notice that increasing the value of λ OH exponentially reduces the number of solutions not-respecting the constraint.Nevertheless, for high values of λ OH , we can notice a reduction of the segmentation accuracy.This reduction is due to the increased energy of the potential barriers in the quantum annealer between two solutions respecting the one-hot encoding constraint.As quantum annealers use tunneling effect to converge to the global optimum, increasing the energy of these barriers reduces the probability of the system to converge to a lower minimum [33,40].Consequently, even if all tested values of λ OH respect the requirements presented in section 3.1, its value has to be judiciously chosen to limit the error rate whilst maximizing the segmentation accuracy.For error labelled pixels, completions methods based on the labels of nearest neighbours can be applied to achieve a satisfying segmentation of the image, without loosing the quantum speedup advantage.
Concerning the annealing time, Fig. 7 also shows a correlation between the error rate and the annealing time, which corroborates the results of [40].On this figure, we have chosen a value of λ OH such that the impact of the annealing time on the number of errors is clearly visible (same for Fig. 6).Optimal annealing time is around 250µs which is usually a good value [40] confirmed by our experiments.Disregarding the error of the quantum hardware, larger instances may be implemented by duplicating some variables with a lot of connections.Even if D-Wave embedding functions already binds logical variables to multiple physical qubits, these generic functions may not provide an optimal embedding.Hence, including the optimal embedding in the QUBO formulation could increase the number of implementable instances, at the cost of additional ancillary variables.Moreover, one major hypothesis of the model is the number of classes in the image.In practice, for a SAR image, there is no guarantee that a vehicle is present on the image.Then, the image should have to be segmented in a single "background" class to maximize the accuracy, whilst the EM algorithm will try to find 3 distinct classes.Hence, it could be accurate to use the EM algorithm to find optimal parameters on reference images (with 3 classes) in order to have an initial guess of the classes parameters, then execute the MRF with fixed parameters and Q = 3 on all the images of the dataset, including images without vehicles.For similar noise parameters, this approach would achieve a satisfying segmentation of all images, including images without vehicle.Another method could be to reduce the number of classes during the EM algorithm process.At a given epoch, if the number of pixels labelled in a certain class is below a predefined limit, we could remove this class and proceed the algorithm for a reduced number of classes.However, this method appears to be highly error sensitive, as a single irrelevant epoch could lead to a definitive removal of a class.

Conclusion
In conclusion, our work proposes a new hybrid quantum-classical method for image segmentation, taking profit of the quantum annealing speedup to perform fast segmentation of SAR images.We also present our first qualitative and quantitative results on the D-Wave quantum annealer and evaluate the scalability of our method regarding the current hardware constraints.Moreover, we propose a discussion on the optimal parameters for the quantum computation, especially on the weighting of constraint cost functions and the annealing time.We conclude that quantum approaches for image segmentation are promising, but improvements of current hardware limitations are required, especially concerning the error rate and the embedding capability.
To complete this paper, further studies could be done on other hardware related parameters, such as optimal embedding.Concerning the hardware, as our formulation is compatible with any quantum annealers, our method could be tested on other quantum computers, as the Pasqal machine [41], which uses spatially arranged cold atoms to perform quantum annealing.Our work can also be extended to other image processing problems useful in the radar domain, as edge detection or denoising.

Figure 1 :
Figure 1: Variation of the accuracy depending on the value of λP .These results are obtained by simulating the multi-class segmentation approach on a (40, 40) image with Q = 4.

Figure 2 :
Figure 2: Qualitative results of the EM algorithm for images with additional Gaussian noise: (a) and (e) corresponds the original image, (b) and (f) corresponds to the original image with the addition of a gaussian noise, (c) and (g) correspond to the segmentation obtained with the simulated annealing approach.(d) and (h) corresponds to the segmentation obtained with the region smoothing algorithm [2]

Figure 3 :
Figure 3: Example of MSTAR SAR image: (1) corresponds to the vehicle, (2) corresponds to the shadow of the vehicle, (3) corresponds to the (noisy) background of the image.

Figure 4 :
Figure 4: Qualitative results of the EM algorithm on M-STAR SAR images.(a/e/i) correspond to the original image.(b/f/j) correspond to the results obtained with our approach.(c/g/k) corresponds to the results obtained with the region smoothing algorithm.(d/h/l) corresponds to the results obtained with the fuzzy clustering algorithm.

Figure 5 :
Figure 5: Qualitative results of the EM algorithm on SSDD dataset images.(a/e/i) correspond to the original image.(b/f/j) correspond to the results obtained with our approach.(c/g/k) corresponds to the results obtained with the region smoothing algorithm.(d/h/l) corresponds to the results obtained with the fuzzy clustering algorithm.

Figure 6 :
Figure 6: Variation of the segmentation accuracy and the percentage of pixels not respecting the one-hot encoding constraint depending on the value of λOH .These results are obtained by executing the multi-class segmentation approach on a (15, 15) image with 2 classes, for 100 samples and an annealing time of 50µs.

Figure 7 :
Figure 7: Variation of the percentage of pixels not respecting the one-hot encoding constraint depending on the annealing time.These results are obtained by executing the multi-class segmentation approach on a (15, 15) image with 2 classes, for 100 samples and λOH = 5.

Table 2 :
Comparison of average computation time and accuracy between simulated annealing and other nonquantum approaches proposed in