Q-Seg: Quantum Annealing-Based Image Segmentation

Introduction

In this post, we’ll explore Q-Seg, a quantum-powered unsupervised image segmentation algorithm designed to solve segmentation challenges in complex datasets like satellite imagery, where traditional methods fall short. This post will delve into both the mathematics behind Q-Seg and a practical implementation using a Jupyter notebook.

Why Quantum for Image Segmentation?

Classical segmentation methods rely on extensive labeled data and can struggle with noisy or ambiguous regions. Q-Seg leverages quantum computing’s ability to explore vast solution spaces rapidly, achieving segmentation without requiring labeled data or extensive preprocessing.

---

The Math Behind Q-Seg

1. Image to Graph Representation

The first step in Q-Seg is representing the image as a graph. As an image is a grid of pixel intensity values, an equivalent graph representation is when:

• Each pixel becomes a node.
• Edges between nodes carry weights based on pixel similarity.

In our method, edge weights, $w(p_i, p_j)$, are calculated with the Gaussian similarity metric:

\[w'(p_i, p_j) = 1 - \exp\left(-\frac{(I(p_i) - I(p_j))^2}{2\sigma^2}\right)\]

This metric yields values between 0 and 1, where higher values indicate greater similarity. These values are then normalized to a range of ([-1, 1]).

\[w(p_i, p_j) = -1 \times \left( \frac{(b-a) \cdot (w'(p_i, p_j) - \min(w))}{\max(w) - \min(w)} + a \right)\]

where $a=-1$ and $b=1$ are the desired bounds. This formulation of the edge weights provides a general similarity score that effectively works, but every specific use case might require a custom edge weight metric.

2. Minimum Cut Formulation as a QUBO Problem

The goal is to divide the graph into distinct regions by cutting the least-similar edges. In graph terms, we’re finding the “minimum cut” that segments the image efficiently. This cut problem is translated into a Quadratic Unconstrained Binary Optimization (QUBO) problem, compatible with quantum annealing:

\[x^* = \arg \min_x \sum_{1 \leq i < j \leq n} x_{v_{i}} (1 - x_{v_{j}}) w(v_{i},v_{j})\]

where $x$ is a binary vector representing which nodes belong to each segmented region.

3. Solving with Quantum Annealing

Quantum annealing, and specifically the D-Wave Advantage, efficiently solves the QUBO formulation through quantum tunneling. Here’s a quick summary of the annealing process:

• Initialization: The system begins in a superposition of all possible states.
• Annealing: Quantum tunneling and entanglement allow the system to “explore” multiple solutions simultaneously.
• Final State: The lowest-energy state corresponds to the optimal segmentation mask.

This quantum approach can quickly navigate the exponentially large solution space that would be computationally intense for classical solvers, especially as image sizes grow.

---

Experimentation and Implementation Guide

Now that we’ve covered the math and the algorithm’s workflow, let’s walk through the implementation in Python. We’ll use a Jupyter notebook to illustrate each step and run an experiment with a real quantum annealer on cloud.

Implementation

This guide provides instructions for setting up and running the Q-Seg project, a quantum-based image segmentation algorithm leveraging D-Wave’s quantum annealer. Follow these steps to get started.

Prerequisites

• Python Version: Has been tested on Python 3.9.11 and above
• Jupyter Notebook: For interactive examples and visualization.

Installation

1. Clone the Repository: Begin by cloning the Q-Seg repository to your local machine. Open a terminal/cmd/powershell and run

   git clone https://github.com/supreethmv/Q-seg.git
   cd Q-seg
   

2. Set Up a Virtual Environment (Optional): Following good practices and to avoid dependency conflicts, create and activate a virtual environment.

   python3 -m venv qseg_env # or `virtualenv qseg_env`
   source qseg_env/bin/activate  # On Windows, use `qseg_env\Scripts\activate`
   

3. Install Dependencies: Use the provided requirements.txt to install all necessary libraries. The process may take a few minutes depending on your internet speed and the packages already available in your system.

   pip install -r requirements.txt
   

Code Walkthrough - Tutorial

After setting up the environment, open the Q-Seg Tutorial Jupyter notebook in the root of the repository directory to start exploring the operational pipeline.

   jupyter notebook tutorial.ipynb

---

Import Dependencies

from warnings import simplefilter
simplefilter(action='ignore', category=FutureWarning)

# Importing functions from the modules in the qseg package
from qseg.graph_utils import image_to_grid_graph, draw, draw_graph_cut_edges
from qseg.dwave_utils import annealer_solver
from qseg.utils import decode_binary_string

# Additional necessary imports
import numpy as np
import matplotlib.pyplot as plt
import networkx as nx
import time

Create a sample Image

height,width = 3,3

image = np.array([
       [0.82,  0.1, 0.99],
       [0.83,  0.2, 0.95],
       [0.1,  0.05, 0.98]
       ])
plt.imshow(image, cmap=plt.cm.gray)

<matplotlib.image.AxesImage at 0x1bc48a1a2d0>

Convert Image to Graph

normalized_nx_elist = image_to_grid_graph(image)  # We are using Guassian similarity metric as the edge weight metric
G = nx.grid_2d_graph(image.shape[0], image.shape[1])
G.add_weighted_edges_from(normalized_nx_elist)
draw(G, image)

Solve using D-Wave annealer

Note: This requires you to have an account in D-Wave Ocean, up on logging in, the D-Wave Leap dashboard is displayed, you can find the private token on the left column under Solver API Token

If you do not wish to run the problem on the remote Quantum Annealer, you can also try out the Gurobi solver which runs locally on your machine. Scroll down for code.

start_time = time.time()
samples_dataframe, execution_info_dict = annealer_solver(G, private_token = 'xxx-xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx')
total_time = time.time()-start_time
execution_info_dict['total_time'] = total_time
execution_info_dict

{'qpu_sampling_time': 182880.0,
 'qpu_anneal_time_per_sample': 20.0,
 'qpu_readout_time_per_sample': 50.86,
 'qpu_access_time': 198804.76,
 'qpu_access_overhead_time': 679.24,
 'qpu_programming_time': 15924.76,
 'qpu_delay_time_per_sample': 20.58,
 'total_post_processing_time': 95.0,
 'post_processing_overhead_time': 95.0,
 'problem_formulation_time': 1.258246898651123,
 'connection_time': 0.18280720710754395,
 'embedding_time': 0.03640317916870117,
 'response_time': 0.023007631301879883,
 'sample_fetch_time': 0.001997709274291992,
 'total_time': 4.44229531288147}

What actually happened?

Figure: Operational pipeline of Q-Seg using the D-Wave quantum annealer.

The above figure gives an overview of all the processes that took place in order to solve the segmentation problem as a QUBO problem on the D-Wave Quantum Annealer.

All samples received

samples_dataframe

	0	1	2	3	4	5	6	7	8	chain_break_fraction	energy	num_occurrences
0	0	1	0	0	1	0	1	1	0	0.0	-4.0870	709
1	1	0	1	1	0	1	0	0	1	0.0	-4.0830	1288
2	1	1	0	1	1	0	1	1	0	0.0	-2.5878	1
3	0	0	1	0	0	1	0	0	1	0.0	-2.5720	1
4	0	1	1	0	1	1	1	1	0	0.0	-1.5098	1

Decode binary solution string to a binary segmentation mask

solution_binary_string = samples_dataframe.iloc[0][:-3]
segmentation_mask = decode_binary_string(solution_binary_string, height, width)
plt.imshow(segmentation_mask, cmap=plt.cm.gray)

<matplotlib.image.AxesImage at 0x1f6b6f54dd0>

cut_edges = [(u, v) for (u, v, d) in G.edges(data=True) if segmentation_mask[u]!=segmentation_mask[v]]
cut_edges

[((0, 0), (0, 1)),
 ((0, 1), (0, 2)),
 ((1, 0), (2, 0)),
 ((1, 0), (1, 1)),
 ((1, 1), (1, 2)),
 ((2, 1), (2, 2))]

draw_graph_cut_edges(G, image, cut_edges)

Solve using Gurobi Solver (Classical State-of-the-art)

import gurobipy as gp
from gurobipy import GRB
from qiskit_optimization.applications import Maxcut # Used just for formulating the Maxcut problem as a QUBO

def gurobi_qubo_solver(G):
  w = -1 * nx.adjacency_matrix(G).todense()
  max_cut = Maxcut(w)
  qp = max_cut.to_quadratic_program()
  linear = qp.objective.linear.coefficients.toarray(order=None, out=None)
  quadratic = qp.objective.quadratic.coefficients.toarray(order=None, out=None)

  linear = {int(idx):-round(value,2) for idx,value in enumerate(linear[0])}
  quadratic = {(int(iy),int(ix)):-quadratic[iy, ix] for iy, ix in np.ndindex(quadratic.shape) if iy<ix and abs(quadratic[iy, ix])!=0}

  qubo_matrix = np.zeros([len(linear),len(linear)])
  for key,value in linear.items():
    qubo_matrix[int(key),int(key)] = value
  for key,value in quadratic.items():
    qubo_matrix[int(key[0]),int(key[1])] = value/2
    qubo_matrix[int(key[1]),int(key[0])] = value/2

  n = qubo_matrix.shape[0]
  model = gp.Model()
  x = model.addVars(n, vtype=GRB.BINARY)
  obj_expr = gp.quicksum(qubo_matrix[i, j] * x[i] * x[j] for i in range(n) for j in range(n))
  model.setObjective(obj_expr)
  model.setParam('OutputFlag', 0)
  model.optimize()

  if model.status == GRB.OPTIMAL:
    solution = [int(x[i].X) for i in range(n)]
    binary_string = ''.join(str(bit) for bit in solution)
    return binary_string, model.objVal
  else:
    return None, None

start_time = time.time()
gurobi_qubo_solution, gurobi_qubo_value = gurobi_qubo_solver(G)
total_time = time.time()-start_time

Set parameter Username
Academic license - for non-commercial use only - expires 2025-02-25

segmentation_mask = decode_binary_string(gurobi_qubo_solution, height, width)
plt.imshow(segmentation_mask, cmap=plt.cm.gray)

<matplotlib.image.AxesImage at 0x1bc48d5a210>

NOTE: Gurobi is an off-the-shelf solver, it is licensed and not open source. The free version allows to run problems upto size 200 only, which is the number of nodes in our case. However, if you are a student or part of an educational institute, you can obtain an academic license for non commercial use for free. Check this tutorial by Gurobi on how to obtain an academic license if you are a student, researcher or faculty.

Here’s a video for a quick overview of the method and results

References

• Paper: IEEE Computer Graphics and Applications
• arXiv: 2311.12912
• Code: Official GitHub Repo