Cannyのエッジ検出から角度を推定するアルゴリズムは以下だと知られている。

対象画像にガウシアンフィルタをかける
縦横のSobelフィルタをかける
着目ピクセルにおける縦/横Sobelの出力を$g_y,g_x$とすると、$\arctan \frac{g_y}{g_x}$が角度値になる

なお、3で出てくる角度は縦線($g_y$出力が0)が角度0に対応することに注意する。

ところで、このなぜこの処理を行うと角度がでるのかを説明している記事が少ない気がするので、この記事で導出および検証してみよう。

導出

一旦、処理対象が画像であることは忘れて直線に沿って非ゼロ、あとはゼロである関数を説明モデルとして定義しよう。

直線角度が$\theta$で表現されるとき、

$$x=y\tan\theta$$

を満たす場合のみ非ゼロとなるような2変数関数を作ればいい。ここで、角度0が縦線に対応するように式を立てている。

$$f(x,y)=\delta(x-y\tan\theta)$$

ここで$\delta$はディラックのデルタ関数であり、以下の特性をもつ。

$$\delta(x)=\begin{cases}
\infty & \text{if $x=0$} \
0 & \text{otherwise}
\end{cases}
$$

$$\int_{x_0}^{x_1}\delta(x-a)f(x)dx=f(a)$$

ただし、$x_0<a<x_1$とする。

いま、これにガウシアンフィルタを畳み込む。

$$g(x,y)=\frac{1}{Z}\int\int f(x-x',y-y')\exp(-x'^2-y'^2)dx'dy'$$

$$g(x,y)\propto \int \exp(-(x-(y-y')\tan\theta)^2-y'^2)dy'$$
$$=\int \exp(-y'^2(1+\tan^2\theta)-2(x-y\tan\theta)y'\tan\theta-(x-y\tan\theta)^2)dy'$$
$$=\int \exp(-(1+\tan^2\theta)(y'+\frac{(x-y\tan\theta)\tan\theta}{1+\tan^2\theta})^2-\frac{(x-y\tan\theta)^2}{1+\tan^2\theta}(1+\tan^2\theta-\tan^2\theta))dy$$

$$=\int \exp(-A(y'-B)^2)dy'\exp(-\cos^2\theta(x-y\tan\theta)^2)$$

$$=W\exp(-U(x-y\tan\theta)^2)$$
ここで、$A,B,W,U$は$x,y$には依存しない($\thetaには依存する$)。

ガウシアンフィルタは分散も指定できるが、適切な$U$に置き換えれば式の形は同じである。

$$g(x,y)=W\exp(-U(x-y\tan\theta)^2)$$

いま、Sobelフィルタの役割をx,yそれぞれの偏微分に対応していると考える。それぞれの偏微分は以下。

$$g_x=\frac{\partial g}{\partial x}=-2UW(x-y\tan\theta)g(x,y)$$
$$g_y=\frac{\partial g}{\partial y}=2UW\tan\theta(x-y\tan\theta)g(x,y)$$

これより、任意の(x,y)において

$$\tan\theta=-\frac{g_y}{g_x}$$

が導ける。

どこがキモか？

Cannyによる角度検出の処理のなかでGaussianフィルタを使う部分があるが、このフィルタがキモだといいたい。
Sobelフィルタが微分に対応しているのは係数の形をみれば自明であるため、重要なのだがつまらない。Gaussianフィルタが以下の特性をもつからこそ、Sobelフィルタの比が$\tan\theta$になるためである。

直線上にGaussianを畳み込むと直線からの距離を変数としたGaussianが得られること
Gaussianを微分すると、距離の方向別偏微分が元の信号にかけ算した結果で出てくること

この特性によって、$\sigma$の値やpixelの位置によらず方向別の偏微分の比が$\tan$に依存するというのが興味深い。

実装と検証

以下のように、画像内に任意角度の一本線を描画して、各角度が上記式によって求まるかを調べてみよう。

import numpy as np
import matplotlib.pyplot as plt

yx=np.indices((32,32)).astype(float)
ct=np.array(yx.shape[1:])/2
yx-=ct[:,None,None]

def makePattern(_deg):
	# ラジアンに変換.
	theta=np.pi*_deg/180
	# 定義式に基づいてdelta(r)のrを生成.
	r=(yx[1]*np.cos(theta)-yx[0]*np.sin(theta))
	# ガウシアンフィルタ後. σは9.
	im=np.exp(-(r**2)/9/9)
	# 縦Sobel.
	ys=(im[2:,:]-im[:-2,:])
	ys=ys[:,:-2]+ys[:,1:-1]*2+ys[:,2:]
	# 横Sobel.
	xs=(im[:,2:]-im[:,:-2])
	xs=xs[:-2]+2*xs[1:-1]+xs[2:]
	# 公式を使ってthetaを逆算. 単位は度に直す.
	angl=-np.arctan2(ys, xs)*180/np.pi
	# 線は-90から90度で一周するので、値域に収まるように180を加減する.
	angl[angl<-90]+=180
	# 90度付近の値の不連続がでにくいように、89.999とする.
	angl[angl>89.999]-=180
	return im, angl

d=[]
resd=[]
for t in np.linspace(-90,90,19):
	#fig=plt.figure()
	#axes=fig.subplots(3,1)
	im,angl=makePattern(t)
	plt.subplot(311).imshow(im)
	plt.subplot(312).imshow(angl,vmin=-90,vmax=90,cmap="jet")
	d.append(t)
	# 理屈上はr=0となる点以外ならどこ使っても同じ値が出るが、線の芯から離れるほど信号が小さくなるので誤差が大きくなる。(1,2)付近の点を角度vs計算結果のプロットの代表点として使用する.
	resd.append(angl[*ct.astype(int)+(2,1)])
	ax=plt.subplot(313).axes
	ax.set_ylim(-90,90)
	ax.set_xlim(-90,90)
	plt.subplot(313).plot(d,resd, marker="o")
	plt.savefig("tan%+03d.png"%t)
	if t<90:
		plt.pause(0.5)
		plt.cla()
		plt.clf()

plt.show()

結果は以下。
それぞれ、一番上の図は指定角度の直線にGaussianを畳み込んだフィルタ画像。中段は、各pixelで角度計算を行った結果である。下段は特定のpixelに着目して、指定角度vs計算結果をプロットしていった結果である(-90から90まで順に描画されていく途中経過)。
中段において中央のpixelなどが0になっているのは、Gaussianの中心点は傾きが0になってしまい角度計算ができないことを意味している。

-40度のとき。

+30度のとき。

+90度のとき。

直線の傾きは+90度とマイナス90度で差がないため、最後は一周して-90が計算結果になっている。
というわけで、導出の通りの結果が得られていることがわかる。

補足。実は一般的に成り立つ話だった

任意の微分可能な関数$h$について

$$g(x,y)=h(x-y\tan\theta)$$

とできるとき、

$$dg/dx=h'(x-y\tan\theta)$$
$$dg/dy=-\tan\theta h'(x-y\tan\theta)$$

$$\tan\theta=-g_y/g_x$$

一瞬で求まってしまった。つまり、やはりGaussianである必要はなく、任意の微分可能関数で定義された畳み込みフィルタを使えばOKということらしい。ややこしく考えすぎただけのようだ。