Hough変換と、局所的な傾きを用いたその拡張

# はじめに
画像から**直線や円を検出**する手法として最もポピュラーなものの一つに**Hough変換**があります。「**ハフ変換**」と読むそうです。

Hough変換は最も**シンプルで理解しやすい図形検出の手法**であり、現在に至るまで**さまざまな改良**が施されています。
初学者向けに、**Hough変換の流れ**と、局所的な傾きを用いた**Hough変換の拡張**について解説します。

# Hough変換の流れ
## 概要
簡単のため、ここからはまず**Hough変換によって直線を検出**する手法について解説します。
画像は以下のように**エッジ化**されているものとします。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/ae9b8ed3-3e54-ba2e-0b0c-7449d2678d13.png" width="35%"><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/69de8c33-a776-5a61-c3d6-8ad2904a4fbe.png" width="35%">
例えば上記の右側ような画像から**直線を検出**することがゴールです。

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/c7f611fb-45b1-f637-bdc4-49b0f4dcc8c4.png" width="35%">

上の画像では、**ノートの端**と**パソコンの端**、**タブレットの端**がそれぞれ正常に検出されていることが確認できます。

## 下準備する
まずは、直線を扱いやすくするために画像を**座標平面上**に載せます。
$x,\ y$の範囲は、ここではどちらも$-1$から$1$までとしています。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/39e0372a-0dee-6b1e-7fc5-958d05b167fe.png" width="35%">

## 投票をとる
次に、**画像のピクセル**を一つ一つ見ていきます。

エッジ上にないピクセルは**スルー**します。
あるピクセル$\mathrm{P}_0$が**エッジ上にあるのを検出**したとしましょう。
以下の図の第一象限にある**黄色い点**が$\mathrm{P}_0$です。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/0eab7e4f-7cfb-b490-1015-2d4ab0e6ca82.png" width="35%">
ピクセル$\mathrm{P}_0$に相当する**座標**を$(x_0,\ y_0)$とします。

ここで、検出する**直線の式**を
$$x\cos\theta+y\sin\theta=r\tag{1}$$
とおきます。$y$軸に平行な直線も統一的に処理するための表し方です。$\tan\theta\neq0$のもとでこれを$y$について解くと
$$y=\frac{1}{\tan\theta}x+\frac{r}{\sin\theta}$$
となるので、$\tan\theta=0$も含めれば$(1)$で**任意の直線の式を網羅**できていることがわかるかと思います。

さて、今**エッジ上の点**$\mathrm{P}_0$を見つけ、そこに相当する**座標**が$(x_0,\ y_0)$で与えられたのでした。
このエッジ上の点を通る**直線$l$の存在を仮定**すると、$l$が$\mathrm{P}_0$を通ることから、$l$の**式**は
$$l:\ x_0\cos\theta+y_0\sin\theta=r\tag{2}$$
と表すことができます。これは、**$\theta$と$r$に関する関係式**と見ることができます。

ここで、以下のような「**投票箱**」を用意します。

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/73a155dd-a9c6-6cf5-af27-6967a4918c27.png" width="35%">

$R$は、設定した座標平面にあわせて適当に設定する**正の実数**です。**検出された直線が座標平面内に収まる**という条件から、三角不等式を用いて定めることができます。今回の場合の値は$2\sqrt{2}$です。

そして、**関係式(2)の条件**を満たす位置に存在する「**投票箱**」に票を入れていきます$^{[1]}$。
なおこのとき、**三角関数の合成**を用いれば(2)式は$\theta r$平面上で**正弦波**を表すことがわかります。

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/714f7d2b-a777-b101-9384-531467ef75b1.png" width="35%">

$\mathrm{P}_0$に相当する座標$(x_0,\ y_0)$を通る直線は、下図の**黄色い直線**のように、**無数に存在**します。ここでの投票は、これら**無数の直線のパラメータについて投票**したことに相当します。

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/81875161-662f-67b4-12cb-a3d61c224fae.png" width="35%">

そして、ここまで紹介した処理を、**エッジ上の各点**について繰り返していきます。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/c9c49cc8-4c87-c925-1bca-056753f764a8.png" width="35%"><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/3abd1e82-b75c-df79-28bd-1f7111e405ef.png" width="35%">

そうすると、**真のパラメータ付近の点**に票が集まります（右上図）$^{[2]}$。

真のパラメータ付近から**ずれた点**にも票が集まりますが、それは**周囲の方向に分散**します。
そのため、**票が集まっている点に着目**すれば、**直線の真のパラメータ$\theta,\ r$が正しく求められる**というわけです。

## 円検出への応用
直線検出の場合と同様に、**Hough変換を用いて円を検出**することもできます。

直線の場合と同様に、**エッジ上の点**$\mathrm{P}_0$が検出され、それに相当する**座標**が$(x_0,\ y_0)$であるとしましょう。

検出する**円の式**を

$$(x-a)^2+(y-b)^2=c^2$$

とすれば、この**エッジ上の点を通る円$C$の式**は

$$C:\ (x_0-a)^2+(y_0-b)^2=c^2$$

すなわち

$$C:\ (a-x_0)^2+(b-y_0)^2-c^2=0\tag{3}$$
となります。今度は$3$次元配列を用意して、(3)式を満たす点$(a,\ b,\ c)$について直線の場合と同様に**投票**を取っていけば、**真のパラメータ付近に票が集中**し、**円を検出**することができます。

なお、この場合パラメータ空間$(a,b,c)$内で(3)式は**円錐面**を表します。
$3$次元に限った話ではありませんが、このパラメータ空間のことを「**Hough空間**」と呼ぶそうです。

## 問題点
Hough変換は**汎用性の高い手法**ですが、**いくつかの問題点**も存在します。

### 計算量が多い
Hough変換は通常**重い処理**です。ピクセルごとに全て投票を取っていると**計算量が膨大**になってしまうため、通常は**投票をとるピクセルを一定間隔で間引く**など何らかの方法で**工夫**を施します。最近では、**投票をとるピクセルを確率的に選ぶ**という手法も盛んに研究されています。

### 隣接する直線の検出におけるジレンマ
画像によっては、画像のノイズや隣接して存在する直線などによって、**複数の投票の山が生成**されることがあります。**近くの山同士を強制的に併合**するようにすれば「同じ直線の上に、**パラメータがわずかに異なる複数の直線が生成される**」という状況を防ぐことができます。しかし同時に「実際には**異なる直線上**にあるのに、**同じ直線として検出されてしまう**」という問題が発生します。なお、投票後の直線の併合には、**クラスタリング**の手法を活かすことができます。

### 次元数が増えると精度が急激に低下する
**Hough変換の致命的な欠点**として、**パラメータ数が多い場合には急激に精度が低下する**という点があります。実際の画像には**ノイズ**も多く含まれるため、**パラメータ付近の広い範囲に点が分散**するようになり、実用上使い物にならなくなってしまうのです。そのための解決策として、次章で紹介する**局所的な傾きを用いる**方法があります。

# 局所的な傾きを用いたHough変換の拡張
ここからは、主に**円を検出**するために**Hough変換を拡張した手法**の解説です。
**独自に考えた手法をベースに解説**しますが、古くから**類似の手法が提案**されているようです$^{[3]}$。

直線の場合はそれほど恩恵がありませんが、わかりやすさのためまずは**直線の場合について解説**します。

## 下準備する
従来のHough変換の際と同様に、まずは**画像に座標平面を設定**します。
その次に、**エッジ上の各ピクセル**に対し、その付近の**局所的な傾き**を求めます。
局所的な傾きは**エッジ検出の過程で副次的に求められることも多い**ですし、**愚直な実装**でも現実的な時間で求まります。
今回の画像の場合、**局所的な傾きを図示**すると以下のようになります$^{[4]}$。

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/1bbb4908-993f-5263-58a7-7da44718791f.png" width="35%">

## 投票をとる
**エッジ上の点**$\mathrm{P}_0$が検出され、相当する**座標**が$(x_0,\ y_0)$であったとします。また、その場所の**局所的な傾きを表す角度**が$\theta_0$であったとします（下図）。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/a00de091-1e2b-24a1-3463-3965ec18d88d.png" width="35%">
（ただし、**$\theta_0$は$-\pi/2<\theta_0\leq\pi/2$の範囲**で、$x$軸から反時計回りの向きに取ります。すなわち上図の場合$\theta_0<0$です）

ここで、**投票箱**を用意します。

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/73a155dd-a9c6-6cf5-af27-6967a4918c27.png
" width="35%">

$\mathrm{P}_0$を通る**直線$l$の存在を仮定**すると、$l$の式は
$$l:\ x_0\cos\theta+y_0\sin\theta=r\tag{4}$$

と表すことができます。これは$\theta$と$r$**に関する関係式**と見ることができますが、$\theta_0$が**既知である**ことに注意すると、これは$(\theta,\ r)$**平面上におけるある一点**$(\theta_0,\ x_0\cos\theta_0+y_0\sin\theta_0)$を表すことがわかります。
したがって、$P_0$に関する**投票を実施**したあとの**投票箱**は以下のようになります。

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/04e92c0a-a853-3d28-dc2c-e5e1e77059c8.png" width="35%">

この投票を**エッジ上の全ての点**について実施すると、**投票実施後の投票箱**は以下のようになります。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/ff428f51-8937-855b-c102-4729d1a46399.png" width="35%">

先に**局所的な情報**を調べておくことで、**無駄な投票を削減**し、**精度を向上**させることができます。

## 円検出への応用
同様に、**円検出へ応用**することも可能です。

**エッジ上の点$\mathrm{P}_0$が検出**され、相当する**座標**が$(x_0,\ y_0)$、**局所的な傾きを表す角度**が$\theta_0$であったとします。

このとき、$\mathrm{P}_0$を通り、$x$軸となす角度が$\theta$で表される**直線の方程式**は
$$(\sin\theta_0)x-(\cos\theta_0)y+y_0\cos\theta_0-x_0\sin\theta_0=0\tag{5}$$
と表されます。
一方、**円の方程式**を
$$C:\ (x-a)^2+(y-b)^2=c^2$$
と表すと、$C$上の点$(x_0,\ y_0)$における**接線の方程式**は
$$(x-a)(x_0-a)+(y-b)(y_0-b)=c^2$$
すなわち
$$(x_0-a)x+(y_0-b)y+a^2+b^2-c^2-ax_0-by_0=0\tag{6}$$
と表されます。**(5)と(6)の係数を比較**して
```math
\begin{cases}
\sin\theta_0(y_0-b)&=-\cos\theta_0(x_0-a)\\
(a^2+b^2-c^2-ax_0-by_0)\sin\theta_0&=(x_0-a)(y_0\cos\theta_0-x_0\sin\theta_0)
\end{cases}
```
を満たす$3$**次元空間内の点**$(a,\ b,\ c)$**に投票**すればよいです。**ゴツそうな見た目**をしていますが、第二式の右辺にある$(y_0\cos\theta_0-x_0\sin\theta_0)$などは**単なる数値**なので、実際はそれほどややこしくありません。

なお、**$3$次元空間内に投票するのはコストがかかる**ため、**次で紹介するような改良**も可能です。

### 円検出における改良
**円の半径と中心の座標を同時に発見するのは大変**なので、**困難を分割**します。すなわち、**先に円の中心**を求め、**後からその円の半径**を求めるという戦略をとります。

**円の中心の座標**を求めます。
エッジ上の点$\mathrm{P}_0$における**局所的な傾きと直交する傾き**を持ち、**点$\mathrm{P}_0$を通る直線**$m$の式は
$$m:\ (\cos\theta_0)x+(\sin\theta_0)y-x_0\cos\theta_0-y_0\sin\theta_0=0$$
と表されます。
ここで、**円の中心を求めるための$2$次元の投票箱**を用意し、$m$上にある投票箱に投票していきます。
**円上の点**から、その点における**傾きと直交する向きに引かれた直線**は**必ず円の中心を通る**ので、円上の点からの投票を重ねると**投票箱は以下**のようになります。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/6fe4fe9d-6771-152a-5baa-ebbe158b6f66.png" width="35%">
ここで投票が重なった点が、**円の中心**です。

その後、**検出された中心とエッジ間の距離を計算**します。ここでは$1$**次元の投票**が実施され、多く投票が集まっているところが**その円の半径**となります。もちろん、同じ点を中心に持ち半径が異なる**複数の円が検出**される可能性もあります。

## おわりに
**画像検出っておもしろい**ですね。**パラメータが多い**場合は、上手くHough空間を定めることで**Hough変換したあとのパラメータをさらにHough変換する**、なんてこともできるんでしょうか。Hough変換は**応用や改良の余地が大きく**、工夫しがいがありました。

### 脚注
[1] ここでは説明のために投票箱を大雑把に区切っています。
[2] 実際はここまで綺麗になりません。また、図においては、黄色い点たちが乗る直線上以外の点からの投票は省略しています。
[3] D.H. Ballard, "Generalizing the Hough Transform to
Detect Arbitrary Shapes", PR, 13-2, 111/122 (1981)
[4] 図ではエッジの周りにも細かい流れがありますが、本来これは除去した方が良い結果が得られます。