はじめに

円を検出するHough変換（ハフ変換）を調べると、後ほど説明する非常に非効率な方法の紹介で止まっていることがほとんどで、説明していても局所勾配を用いた効率化に言及している程度です。

ここでは、位相符号化によるHough空間の２次元化とフーリエ変換の組み合わせによる特徴点の数に依らない計算量の円検出Hough変換を紹介したいと思います。
特にフーリエ変換を用いる方法は他では（何故か）見たことがないので、是非読んでいってください。

Hough変換

Hough変換を知らない方がこの記事にたどり着くことは少ないと思いますが、まず初めに軽くおさらいしておきます。

Hough変換（ハフ変換）とは、画像から直線や円を検出するために用いられる変換で、シンプルでありながら思わず感心してしまうような鮮やかな手法です。

例えば直線に関して例を示すと、
ある点$(x_i,y_i)$を通る直線は無限に存在しますが、それらの直線は原点からの距離とその直線の傾きの角度を用いて
$$r=x_i\sin\theta+y_i\cos\theta$$
と表せます。これは$r-\theta$平面上の正弦波です。
同じ直線上の点であれば、その点を通る直線を表す正弦波は同じ$(r,\theta)$の点を通るはずですので、交点から直線を一意に決定することができるという仕組みです。

図解すると、
次のように直線上に並んでいるような複数の点があるとします。

この点それぞれに上式を適用して$r-\theta$平面上の正弦波として描くと次のようになり、

点の集中度合いを定量化するために2次元ヒストグラムにすると次のようになります。

点の集中度合い上位3つの点に対応するパラメータで直線を描くと以下のように直線が検出できていることが分かります。

ヒストグラムの解像度を上げればパラメータの精度は上がりますが、ノイズに弱くなる点と計算量が増えてしまう点に注意が必要です。

もともとの点が打たれている$x-y$平面から、各点が知りたい図形を表すパラメータ空間（直線の場合は一般的に$r-\theta$平面）に変換することをHough変換といい、このパラメータ空間をHough空間と呼ぶこともあるそうです。

画像から直線などを検出したい場合は、まずエッジ検出などで特徴点を見つけ出して、それら特徴点に対してHough変換を適用します。

円を検出するHough変換

中心座標$(x,y)$と半径$r$の３パラメータが決定できれば円を一意に表すことができます。
ですので、円を検出する愚直なHough変換は$xy$平面→$xyr$空間への変換になります。
ある特徴点を円周上に含む円の中心の候補は、$xyr$空間では特徴点を頂点とした$r$方向の円錐面上となりますので、各特徴点から伸びる円錐面の交点が見出だせれば円を検出することができます。

上の図を見ると円のパラメータは簡単に求まるように思えますが、実際には複数の円やノイズ（円以外の図形を含む）が存在する中から、円錐面が集中する点を発掘しなければいけないため、イメージ的には３次元ヒストグラムのような状態になり非常に重い計算・メモリ消費になる上、パラメータ空間が広い分、精度良く求めるにはより多くの特徴点を必要とします。

ですので計算を削減する工夫として、

ローカルな法線方向を求めて探索方向を絞る
検出対象の半径を限定する

などの方法が取られているようです。

位相符号化用いた円検出Hough変換

実際にこのような非効率に思える方法で計算しているのか気になったため調べていたところ、MATLABの円検出関数の解説で『位相符号化』という気になる文字列を発見しました。
記事内には詳細は記載されていませんでしたが、参照論文¹を見ると具体的な計算が説明されていました。

上記の愚直な方法では、半径を特定するために$r$軸を追加していましたが、位相符号化を使用する方法では半径の表現として位相を使用します。

特徴点からの距離を$R$、検出対象の半径の最大・最小を$R_{max}, R_{min}$とすると、

$$\phi=2\pi\left(\frac{R-R_{min}}{R_{max}-R_{min}}\right)$$

のように位相と半径を対応付けるわけです。²

位相符号化を使用すると何が嬉しいかというと、

探索空間が$x,y$の2次元で済む
同じ円周上以外からの位相は互いに打ち消し合う作用が働く

という非常に大きな利点があります。
計算量・メモリ使用量を抑えられるだけでなく、ノイズにも強くなるのですから嬉しいこと尽くしですね。

フーリエ変換の導入

さて、位相符号化の導入により探索空間が$x,y$のみとなりました。
さらに幸運なことに、円のHough変換は直線のHough変換とは異なり、画像と同じ座標系です。

ここで一旦位相のことは忘れて、直線のHough変換で正弦波を描いた操作が、探索空間を$x,y$に限定した円ではどのような操作に対応するのか考えてみましょう。

直線のときの正弦波の各点は、特徴点を通り得る直線の$(r,\theta)$の組み合わせを表していたのでした。
ですので、円の場合は特徴点を円周上に含み得る円の中心座標$(x,y)$の組み合わせを図示すべきだということが分かります。
探索範囲を$R_{min}\lt R\lt R_{max}$に限定していますので、これは円環になります。

これを全ての特徴点に対して行うと、円の中心を特定することが可能です。

今行った操作を一度立ち止まって整理すると、

全ての特徴点に対して、その特徴点を中心とする$R_{min}\lt R\lt R_{max}$の円環を描画しました。

この操作、何かを思い出しませんか？少し遠いでしょうか？

この操作は畳み込みで置き換えられそうですね。
そして、
畳み込みはフーリエ変換の積（のフーリエ逆変換）と等価なのでした。

フーリエ変換と位相符号化による円検出Hough変換

いよいよフーリエ変換と位相符号化を組み合わせます。

特徴点を検出した画像に円環を畳み込むために、まず双方にフーリエ変換を施します。

\begin{align}
F_{image}(u,v)&=\int{f_{image}(x,y)e^{-2\pi i(ux+vy)}dxdy}\\
F_{annulus}(u,v)&=\int{f_{annulus}(x,y)e^{-2\pi i(ux+vy)}dxdy}
\end{align}

これらをこのまま掛け合わせてフーリエ逆変換すれば前節のようなシンプルな畳み込みの結果が得られますが、
ここで円環のフーリエ変換に位相符号化の処理を加えます。

\begin{align}
\hat{F}_{annulus}(u,v)&=\int{f_{annulus}(x,y)e^{i\phi(x,y)}e^{-2\pi i(ux+vy)}dxdy}\\
&=\int{f_{annulus}(x,y)e^{2\pi i\left(\frac{R-R_{min}}{R_{max}-R_{min}}\right)}e^{-2\pi i(ux+vy)}dxdy}\\
&=\int{f_{annulus}(x,y)e^{2\pi i\left(\frac{\sqrt{x^2+y^2}-R_{min}}{R_{max}-R_{min}}\right)}e^{-2\pi i(ux+vy)}dxdy}
\end{align}

これだけです。
位相符号化とフーリエ変換まで完了したので、これらを掛け合わせてフーリエ逆変換すればHough変換完了です。

\begin{align}
f_{Hough}(x,y)&=\int{F_{image}(u,v)\hat{F}_{annulus}(u,v)e^{2\pi i(ux+vy)}dudv}
\end{align}

実装

Pythonで実装して、実際にどのように機能するのか確認してみます。
数式を見るとなんだか難しそうに思えるかもしれませんが、フーリエ変換も逆変換もnumpyがやってくれるので（自前で実装してもそれほど難しくはありませんが）、自分でやることと言えば円環のマスクを作成することと位相符号化用の$\exp\left[{2\pi i\left(\frac{\sqrt{x^2+y^2}-R_{min}}{R_{max}-R_{min}}\right)}\right]$を用意して掛け合わせておくことくらいです。

今回こちらの画像を使用します。

出典：MathWorks

まず画像から特徴点を抽出します。
今回はOpenCVを使用してCanny法で済ませます。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import cv2

im = plt.imread("画像のパス")
h, w, *_ = im.shape
edge = cv2.Canny(im, 100, 300)
plt.imshow(edge)

あまり綺麗すぎても位相符号化なしとの比較がし辛いのでこの程度にしておきます。

次に円環マスクを作ります。
計算しやすいように一旦原点を中央に持って来ています。

rmin = 22
rmax = 32
h, w = edge.shape
r = ((np.indices((h, w)) - [[[h/2]], [[w/2]]])**2).sum(0) ** .5
mask = (r>=rmin) & (r<=rmax)

plt.imshow(mask)
plt.show()

これを位相符号化せずにそのまま畳み込んでみましょう。

mask_fft = np.fft.fft2(np.fft.fftshift(mask))
edge_fft = np.fft.fft2(edge)
ifft = np.fft.ifft2(mask_fft*edge_fft)
hough = abs(ifft)
plt.imshow(hough, cmap='inferno')

半径の範囲を適切に設定しているのでこれでも中心位置は分かりますね。

続いて位相を符号化します。

mask = mask * np.exp(2j*np.pi*(r-rmin)/(rmax-rmin))

まあ、これだけなんですけどね。

畳み込みの結果を見てみましょう。

mask_fft = np.fft.fft2(np.fft.fftshift(mask))
edge_fft = np.fft.fft2(edge)
ifft = np.fft.ifft2(mask_fft*edge_fft)
hough = abs(ifft)
plt.imshow(hough, cmap='inferno')

先ほどと比べるとSN比が圧倒的ですね。

位相符号化では位相から半径も知ることができるのでした。確認してみましょう。

phase = np.angle(ifft)
phase = np.where(phase>0, phase, 2*np.pi+phase) #-pi~pi -> 0~2pi
r = phase / (2*np.pi) * (rmax-rmin) + rmin
plt.imshow(r*(hough>hough.mean()*7), cmap='inferno', vmin=rmin-5)
plt.colorbar()

少し分かりづらいので、近接点をまとめる処理を足して、実際に円として描画してみましょう。

p = np.array(np.where((hough>(hough.mean()*9))))
l = []
for a, b in np.array(np.where(((p.T[..., None] - p)**2).sum(1)<15)).T:
    for j in l:
        if a in j or b in j:
            j.add(a)
            j.add(b)
            break
    else:
        l.append({a, b})
a = []
for i in l:
    idx = p[:, list(i)]
    a.append(idx[:, acc[idx[0], idx[1]].argmax()])

ax = plt.axes()
ax.imshow(im, cmap='gray')
for y, x in a:
    c = patches.Circle(xy=(x, y), radius=r[y, x], fc='#0000', ec='r', lw=2)
    ax.add_patch(c)

正しい径が得られていそうですね。

これで終了なのですが、最後の円描画時のデータクレンジングはさておき、ハフ変換を行っている部分には一切ループ処理が出現していないことにお気づきでしょうか。
フーリエ変換を用いたHough変換では、特徴点が1個であろうと、数万個であろうと計算量は変化しないのです（画像サイズにのみ依存）。

おわりに

円を検出するHough変換とフーリエ変換はこんなにも親和性が高いのに、誰も触れていない（観測範囲内）ことに疑問を感じて記事に起こしてみました。
OpenCVなどの便利なライブラリが揃っているため、ご自身で実装する機会もあまり多くはないと思いますが、どこかで誰かの参考になれば嬉しいです。

T.J Atherton, D.J. Kerbyson. "Size invariant circle detection." Image and Vision Computing. Volume 17, Number 11, 1999, pp. 795-803. ↩
論文中では線形以外の対応付けも提案されています。 ↩