量子コンピュータで画像のエッジ検出

TL;DR

画像のエッジ検出は、オンラインミーティングでの背景透過や画像認識など、さまざまな応用分野で重要な役割を果たしています。
この画像処理に対して、量子コンピュータは、従来のコンピュータに比べて指数的な優位性を持つため、非常に高効率な処理が可能です。例えば、1000 x 1000ピクセル、64ビットの画像を通常の方法で処理する場合、100万ビット（100万次元の実数ベクトル）が必要ですが、量子コンピュータではこれをわずか20量子ビットで表現できます。
本記事では、このような指数優位性をもつ量子コンピュータにおけるエッジ検出の理論と実装例を紹介します。

対象読者

対象読者としては量子力学にある程度理解のある方になると思います。
テンソル積や観測確率、アダマールゲートなどは解説しません。

理論

古典画像からQPIE形式へ

古典的な画像として2x2ピクセルの次の画像を考えます。画像は白黒で各ビットの色が$I_0$などと書かれているとします。

この画像のビットに対してまず色を次のように規格化し、

\begin{equation}
c_i = \dfrac{I_i}{\sqrt{\sum_j I_j^2}} \tag{1}
\end{equation}

左上から$\ket{00}$、$\ket{01}$、$\ket{10}$、$\ket{11}$と割り当て、

\ket{\mathrm{Img}} = c_0 \ket{00} + c_1 \ket{01} + c_2 \ket{10} + c_3 \ket{11}

と書く状態の線形結合で書くことを考えます。
そうすると、例えば$I_0$の色は$\ket{00}$を観測する確率$|c_0|^2$に対応し、各ピクセルの色が対応する状態の観測確率で表すことができます。
このような各ピクセルをある状態の観測確率で表す形式をQPIE (Quantum Probability Image Encoding)形式といいます。¹

一般の画像を表すためには、Nピクセルの画像の場合には$log_2 N$のqubitが必要となります。(例えば、8x8ピクセルの場合には$log_2 64 = 6$qubitが必要）

QPIE 形式におけるエッジ検出

上記と同様にまず、2 x 2ピクセルの画像で考えます。
$^t[c_0, c_1, c_2, c_3]$に対して次の行列

\begin{align}
I \bigotimes H_0 &=  \begin{bmatrix} 1 & 0 \\0 & 1\end{bmatrix} \bigotimes \dfrac{1}{\sqrt{2}} \begin{bmatrix} 1 & 1 \\ 1 & -1 \end{bmatrix} \\
&= \dfrac{1}{\sqrt{2}} \begin{bmatrix} 1 & 1 & 0 & 0 \\ 1 & -1 & 0 & 0 \\0 & 0 & 1 & 1 \\ 0 & 0 & 1 & -1 \end{bmatrix}
\end{align}

をかけると、

\dfrac{1}{\sqrt{2}} \begin{bmatrix} 1 & 1 & 0 & 0 \\ 1 & -1 & 0 & 0 \\0 & 0 & 1 & 1 \\ 0 & 0 & 1 & -1 \end{bmatrix} \cdot \begin{bmatrix}c_0 \\ c_1 \\ c_2 \\ c_3 \end{bmatrix} = \dfrac{1}{\sqrt{2}} \begin{bmatrix}c_0+c_1 \\ c_0-c_1 \\ c_2+c_3 \\ c_2-c_3 \end{bmatrix}

となります。
得られたベクトルをよく見ると奇数番目に$c_0-c_1$などの色の差分が現れています。

これでは、$c_0$と$c_1$、$c_2$と$c_3$の差分はわかりますが、同様の横方向の$c_1$と$c_2$²や縦方向の差分はわかりません。
まず、横方向ですが、$^t[c_0, c_1, c_2, c_3]$を$^t[c_1, c_2, c_3, c_0]$とすることで得られます。
次に縦方向は画像の転置、すなわち

\ket{\mathrm{Img}} = c_0 \ket{00} + c_2 \ket{01} + c_1 \ket{10} + c_3 \ket{11}

となるような場合を考えると、$^t[c_0, c_1, c_2, c_3]$が$^t[c_0, c_2, c_1, c_3]$となり、$c_0$と$c_1$、$c_1$と$c_3$の差分が得られることになります。

一般のn qubitで表されるような場合にはそれぞれ係数を$^t[c_0, c_1, \cdots, c_{N-1}]$($N=2^n$)としてこのベクトルに

I_{2^{n-1}} \bigotimes H_0 =
\begin{bmatrix}
1 &  1 & 0 &  0 & \cdots & 0 &  0 \\
1 & -1 & 0 &  0 & \cdots & 0 &  0 \\
0 &  0 & 1 &  1 & \cdots & 0 &  0 \\
0 &  0 & 1 & -1 & \cdots & 0 &  0 \\
\vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\
0 &  0 & 0 &  0 & \cdots & 1 &  1 \\
0 &  0 & 0 &  0 & \cdots & 1 & -1 \\
\end{bmatrix}

を作用させることで、

\dfrac{1}{\sqrt{2}}
\begin{bmatrix}
c_0 + c_1 \\
c_0 - c_1 \\
\vdots \\
c_{N-1} + c_{N-1} \\
c_{N-1} - c_{N-1} 
\end{bmatrix}

として差分が得られます。

QHED (Quantum Hadamard Edge Detection)

上記で差分は検出することができましたが、1方向に2回、合計4回に分けて計算する必要がありました。
これを1方向に1回ずつ、合計2回に減らす方法がQHEDと呼ばれる手法です。

まず$\ket{\mathrm{Img}}$に加えて、$H_0\bigotimes \ket{0} = \frac{1}{\sqrt{2}} \ket{0}+\ket{1}$という状態をもつことを考えるとその状態ベクトルは

\ket{\mathrm{Img}}\bigotimes \dfrac{1}{\sqrt{2}} \\(\ket{0}+\ket{1}\\) = \dfrac{1}{\sqrt{2}} \begin{bmatrix} c_0 \\ c_0 \\ c_1 \\ c_2 \\ \vdots \\ c_{N-1} \\ c_{N-1}\end{bmatrix}

と表されます。
この状態に対して、

D_{2^{n+1}}=\begin{bmatrix}
0 & 1 & 0 & \cdots & 0 \\
0 & 0 & 1 & \cdots & 0 \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & 0 & 0 & \cdots & 1 \\
1 & 0 & 0 & \cdots & 0 
\end{bmatrix}

を作用させると、$D_{2^{n+1}} \cdot \ ^t[c_0, c_0, \cdots, c_{N-1}, c_{N-1}] = \ ^t[c_0, c_1, c_1, \cdots, c_{N-1}, c_{N-1}, c_0]$となります。
これに対して、$I_{2^n} \bigotimes H_0$を作用させることで

\dfrac{1}{\sqrt{2}}\begin{bmatrix}
c_0 + c_1 \\
c_0 - c_1 \\
c_1 + c_2 \\
c_1 - c_2 \\
\vdots \\
c_{N-2} + c_{N-1} \\
c_{N-2} - c_{N-1} \\
c_{N-1} + c_{0} \\
c_{N-1} - c_{0}
\end{bmatrix}

となり奇数項で$c_0$から$c_{N-1}$までの各項の差分が奇数行で得ることができました。
縦方向については同様に、画像の転置をとってから同様の処理をすることで求めることができ、最後にそれを統合することで画像のエッジ検出ができます。

実装

実行環境は以下になります。

OS: MacOS 14.2.1
python: 3.12.3
qiskit: 1.1.0

準備

まず必要なパッケージをインストールします。

from qiskit import QuantumCircuit
from qiskit_aer import AerSimulator
from qiskit.visualization import *

import numpy as np
import matplotlib.pyplot as plt

画像のプロットに必要な関数を定義しておきます。

def plot_image(img, title: str):
    plt.title(title)
    plt.xticks(range(img.shape[0]))
    plt.yticks(range(img.shape[1]))
    plt.imshow(img, extent=[0, img.shape[0], img.shape[1], 0], cmap='viridis')
    plt.show()

今回用いる画像を定義します。

image = np.array([[0, 0, 0, 0, 0, 0, 0, 0],
                  [0, 1, 1, 1, 1, 1, 0, 0],
                  [0, 1, 1, 1, 1, 1, 1, 0],
                  [0, 1, 1, 1, 1, 1, 1, 0],
                  [0, 1, 1, 1, 1, 1, 1, 0],
                  [0, 0, 0, 1, 1, 1, 1, 0],
                  [0, 0, 0, 1, 1, 1, 1, 0],
                  [0, 0, 0, 0, 0, 0, 0, 0]])
    
plot_image(image, 'Original Image')

正規化

$(1)$式をもとにピクセルの色を正規化し係数$c$のベクトルを作成します。
画像の転置をとり、横方向に正規化したベクトルと、縦方向に正規化したベクトルとを用意します。

def amplitude_encode(img_data):
    rms = np.sqrt(np.sum(np.sum(img_data**2, axis=0)))
    img_norm = []

    for arr in img_data:
        for ele in arr:
            img_norm.append(ele/rms)
    
    return np.array(img_norm)

image_norm_h = amplitude_encode(image)
image_norm_v = amplitude_encode(image.T)

量子回路

今回必要なqubitは$8\times8=64=2^6$ピクセルであるから6 qubitです。
加えてQHEDで1qubit必要なので、合計7ビット必要です。
また、$D_{2^{n+1}}$を定義しておきます。

data_qb = 6
anc_qb = 1
total_qb = data_qb + anc_qb

D2n_1 = np.roll(np.identity(2**total_qb), 1, axis=1)
# [[0 1 0 ... 0 0 0]
# [0 0 1 ... 0 0 0]
# [0 0 0 ... 0 0 0]
# ...
# [0 0 0 ... 0 1 0]
# [0 0 0 ... 0 0 1]
# [1 0 0 ... 0 0 0]]

実際に量子回路を作成します。

qc_h = QuantumCircuit(total_qb)
qc_h.initialize(image_norm_h, range(1, total_qb))
qc_h.h(0)
qc_h.unitary(D2n_1, range(total_qb))
qc_h.h(0)
qc_h.save_statevector()
display(qc_h.draw('mpl', fold=-1))

qc_v = QuantumCircuit(total_qb)
qc_v.initialize(image_norm_v, range(1, total_qb))
qc_v.h(0)
qc_v.unitary(D2n_1, range(total_qb))
qc_v.h(0)
qc_v.save_statevector()
display(qc_v.draw('mpl', fold=-1))

circ_list = [qc_h, qc_v]

シミュレーション

この回路は実機では実行することができません。
AerSimulatorを用いてシミュレーションを実行し、最終的な状態ベクトルを取得します。

backend = AerSimulator()
backend.set_options(method='statevector')

result = backend.run(circ_list).result()

sv_h = result.get_statevector(qc_h)
sv_v = result.get_statevector(qc_v)

適当な閾値を設定し、奇数部分のみデータをとります。

threshold = lambda amp: (amp > 1e-15 or amp < -1e-15)

edge_scan_h = np.abs(np.array([1 if threshold(sv_h[2*i+1].real) else 0 for i in range(2**data_qb)])).reshape(8, 8)
edge_scan_v = np.abs(np.array([1 if threshold(sv_v[2*i+1].real) else 0 for i in range(2**data_qb)])).reshape(8, 8).T

plot_image(edge_scan_h, 'Horizontal scan output')
plot_image(edge_scan_v, 'Vertical scan output')

最後に縦と横の場合を合わせることでエッジ検出された画像を得ます。

edge_scan_sim = edge_scan_h | edge_scan_v


plot_image(image, 'Original image')
plot_image(edge_scan_sim, 'Edge Detected image')

参考

1. Xi-Wei Yao et al. “Quantum Image Processing and Its Application to Edge Detection: Theory and Experiment”. In: Physical Review X (2017). arXiv: 1801.01465. ↩

2. 今回の2x2ピクセルの場合は$c_1$と$c_2$は斜めの点であり必要ではないとも考えられますが、一般のn x nの場合を考えれば必要となります。 ↩