【数値情報処理b】第2回 輪郭抽出

画像処理で輪郭を抽出することを輪郭抽出（エッジ抽出）という。

画像で輪郭というと物体と物体の境目、あるいは物体と背景の境目である。 画像の濃度や色に急な変化があるところが輪郭に見えることを利用して、輪郭を抽出する。

以下では簡単のためにグレースケール画像で説明する。

濃淡の変化

GIMP による輪郭抽出

GIMP で輪郭抽出を行うには、メニューの [フィルター]→[輪郭抽出]→[輪郭] を選択する。

輪郭抽出にはいくつかのアルゴリズムがあり、GIMP では以下が選べるようになっている。

畳み込み演算

GIMP で畳み込み演算（後述）を行うには、メニューの [フィルター]→[汎用]→[コンボリューション行列] を選択する。

畳み込み演算（コンボリューション）

輪郭抽出などの画像処理の演算には畳み込み演算（コンボリューション）が基本となる。 畳み込み演算とは、ある画素をその周辺の画素の画素値（情報）を利用して処理する演算方法である。

畳み込み演算では、その画素と周辺の画素の画素値に何らかの重みをかけて足し合わせることで変換を行う。 この重みの行列をフィルタ1)と呼ぶ。

例えば、フィルタ $f$ を $x$ 方向が 3個、$y$ 方向が 3個の要素を持つ行列とする。 座標 $(x,y)$ の画素において、画素値 $g$ とフィルタ $f$ の畳み込み演算の結果 $(g*f)(x,y)$ は以下で定義される。

\[ (g*f)(x,y)=\sum_{i=1}^{3}\sum_{j=1}^{3}g(x+i-2,y+j-2)f(i,j) \] \begin{align*} \small &=\ \ g(x-1,y-1)f(1,1)&&+g(x,y-1)f(2,1)&&+g(x+1,y-1)f(3,1)\\ &\ \ \ +g(x-1,y\ \ \ \ \ \ \ )f(1,2)&&+g(x,y\ \ \ \ \ \ \ )f(2,2)&&+g(x+1,y\ \ \ \ \ \ \ )f(3,2)\\ &\ \ \ +g(x-1,y+1)f(1,3)&&+g(x,y+1)f(2,3)&&+g(x+1,y+1)f(3,3) \end{align*} ここで $f(i,j)$ はフィルタ $f$ の要素で、$g(x,y)$ は座標 $(x,y)$ の画素の画素値である。 $(g*f)(x,y)$ が座標 $(x,y)$ の新しい画素値となる。

これは数学における行列の積の定義と異なるので注意してほしい。

畳み込み演算

correlate

imager では畳み込み演算は correlate() で行う。 

> correlate(cimg変数, フィルタ)

同じような関数に convolve() もあるが、上述の説明とフィルタの向きが逆になる。

フィルタは as.cimg() で作成できる。

as.cimg

as.cimg() はベクトルや行列を cimg 型に変換するための関数である。

例えば、$3\times 3$ のフィルタ \[ \left( \begin{array}{ccc} f(1,1)&f(2,1)&f(3,1)\\ f(1,2)&f(2,2)&f(3,2)\\ f(1,3)&f(2,3)&f(3,3) \end{array} \right) = \left( \begin{array}{ccc} 0.1&0.2&0.3\\ 0.4&0.5&0.6\\ 0.7&0.8&0.9 \end{array} \right) \] を作成したい場合は以下を実行する（imager（R のパッケージ） を事前に読み込んでおく必要がある）。 

> f <- as.cimg(c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9), x = 3, y = 3)

ちゃんと作成できたかは以下で確認できる。 

> t(f[,,1,1])
     [,1] [,2] [,3]
[1,]  0.1  0.2  0.3
[2,]  0.4  0.5  0.6
[3,]  0.7  0.8  0.9

t() は行列の行と列の成分を入れ替える転置を行う関数である。 f をそのまま表示すると行列としての表示になり、縦方向が $x$ 方向、横方向が $y$ 方向になってしまうため、転置をして見やすくしている。

オプションとして cc を指定することでカラーチャンネルを加えることもできる。

作成したフィルタ f を使って次の画像データ g と畳み込み演算してみよう。 \[ \left( \begin{array}{ccc} g(1,1)&g(2,1)&g(3,1)\\ g(1,2)&g(2,2)&g(3,2)\\ g(1,3)&g(2,3)&g(3,3) \end{array} \right) = \left( \begin{array}{ccc} 0.1&0.1&0.1\\ 0.2&0.2&0.2\\ 0.3&0.3&0.3 \end{array} \right) \] 

> g <- as.cimg(c(0.1, 0.1, 0.1, 0.2, 0.2, 0.2, 0.3, 0.3, 0.3), x = 3, y = 3)

この g は以下の画素値を持つデータである。 

> t(g[,,1,1])
     [,1] [,2] [,3]
[1,]  0.1  0.1  0.1
[2,]  0.2  0.2  0.2
[3,]  0.3  0.3  0.3

g と f との畳み込み演算を行うと、座標 $(x,y)=(2,2)$ の演算結果 $(g*f)(x,y)$ は以下になる。 

> correlate(g, f)[2, 2]
[1] 1.08

\[ (g*f)(2,2)=\sum_{i=1}^{3}\sum_{j=1}^{3}g(2+i-2,2+j-2)f(i,j) \] \begin{align*} \small &=\ \ g(1,1)f(1,1)&&+g(2,1)f(2,1)&&+g(3,1)f(3,1)\\ &\ \ \ +g(1,2)f(1,2)&&+g(2,2)f(2,2)&&+g(3,2)f(3,2)\\ &\ \ \ +g(1,3)f(1,3)&&+g(2,3)f(2,3)&&+g(3,3)f(3,3)\\ \\ &=\ \ 0.1\times 0.1&&+0.1\times 0.2&&+0.1\times 0.3\\ &\ \ \ +0.2\times 0.4&&+0.2\times 0.5&&+0.2\times 0.6\\ &\ \ \ +0.3\times 0.7&&+0.3\times 0.8&&+0.3\times 0.9\\ \\ &=1.08 \end{align*}

小テスト①

Moodle Server（非公式）で第2回の小テスト①を受験しなさい。

一階微分による輪郭抽出

簡単のためにグレー画像を考える。 黒い背景に白い円が描かれているとすると、我々が輪郭として認識できるところは円の端である。 そこは白黒の濃度が急激に変化するところになる。

画像のある $y$ における $x$ 方向の画素値を $g(x)$ とする。 濃度が急に変化するところは隣り合った画素の画素値が大きく変化するところである。 そこは画素値の微分（導関数）$\frac{dg(x)}{dx}$ が正で大きくなる山や負で小さくなる谷のところだと考えられる。 そこで画素値の微分の値を求めて山や谷になっているところを輪郭として抽出する。 一度だけ微分を行うことを一階微分 $\frac{dg(x)}{dx}$ という。 二度の微分を行うことを二階微分 $\frac{d^{2}g(x)}{dx^{2}}$ という（後述）。

画素値g(x)の一階微分と二回微分

$g(x,y)$ の微分には、$x$ 方向と $y$ 方向の二つがある。 $x$ 方向の微分をするときには $y$ は固定し、$y$ 方向の微分をするときには $x$ を固定する。 この微分を偏微分という。 偏微分の記号には $d$ の代わりに $\partial$（デル）を使う。

\begin{eqnarray*} x方向の偏微分&：&\frac{\partial g(x,y)}{\partial x}\\ y方向の偏微分&：&\frac{\partial g(x,y)}{\partial y} \end{eqnarray*}

また、$x$ 方向と $y$ 方向の偏微分を成分としたベクトルをグラディエントという。 \[ \nabla g(x,y) = \left( \begin{array}{c} \frac{\partial g(x,y)}{\partial x}\\ \frac{\partial g(x,y)}{\partial y} \end{array} \right) \]

ここで $\nabla$ はベクトルの微分演算子でナブラという。 \[ \nabla = \left( \begin{array}{c} \frac{\partial}{\partial x}\\ \frac{\partial}{\partial y} \end{array} \right) \]

輪郭抽出においては、グラディエントのベクトルの大きさ（ノルム）は輪郭の強さ（強度）を表す。 \[ |\nabla g(x,y)| = \sqrt{\left(\frac{\partial g(x,y)}{\partial x}\right)^{2}+\left(\frac{\partial g(x,y)}{\partial y}\right)^{2}} \]

前進差分

偏微分の定義は \begin{eqnarray*} \frac{\partial g(x,y)}{\partial x}&=&\lim_{h\rightarrow 0}\frac{g(x+h,y)-g(x,y)}{h}\\ \frac{\partial g(x,y)}{\partial y}&=&\lim_{h\rightarrow 0}\frac{g(x,y+h)-g(x,y)}{h} \end{eqnarray*} である。 画像はデジタルなので $x, y$ は 1 ずつしか変化せず、連続的に $h\rightarrow 0$ の極限がとれない。 そこで $h=1$ として偏微分を前進差分で近似する。 \begin{eqnarray*} \frac{\partial g(x,y)}{\partial x}\rightarrow g_{x}(x,y)&=&g(x+1,y)-g(x,y)\\ \frac{\partial g(x,y)}{\partial y}\rightarrow g_{y}(x,y)&=&g(x,y+1)-g(x,y) \end{eqnarray*}

これを畳み込み演算で表すと、$x$ 方向と $y$ 方向はそれぞれ以下のフィルタになる。 \[ F_{fx}= \left( \begin{array}{ccc} f(1,1)&f(2,1)&f(3,1)\\ f(1,2)&f(2,2)&f(3,2)\\ f(1,3)&f(2,3)&f(3,3) \end{array} \right) = \left( \begin{array}{ccc} 0&0&0\\ 0&-1&1\\ 0&0&0 \end{array} \right) \] \[ F_{fy}= \left( \begin{array}{ccc} f(1,1)&f(2,1)&f(3,1)\\ f(1,2)&f(2,2)&f(3,2)\\ f(1,3)&f(2,3)&f(3,3) \end{array} \right) = \left( \begin{array}{ccc} 0&0&0\\ 0&-1&0\\ 0&1&0 \end{array} \right) \]

このフィルタで畳み込み演算を行うとそれぞれ以下のようになり、前進差分になっていることが分かる。 \begin{eqnarray*} F_{fx}:(g*f)(x,y)&=&g(x+1,y)-g(x,y)\\ F_{fy}:(g*f)(x,y)&=&g(x,y+1)-g(x,y) \end{eqnarray*}

後退差分

偏微分は以下のように定義してもよい。 \begin{eqnarray*} \frac{\partial g(x,y)}{\partial x}&=&\lim_{h\rightarrow 0}\frac{g(x,y)-g(x-h,y)}{h}\\ \frac{\partial g(x,y)}{\partial y}&=&\lim_{h\rightarrow 0}\frac{g(x,y)-g(x,y-h)}{h} \end{eqnarray*} $h\rightarrow 0$ の極限では先ほどの定義と同等だが、、 $h=1$ として偏微分を後退差分で近似すると前進差分とは少し異なるものになる。 \begin{eqnarray*} \frac{\partial g(x,y)}{\partial x}\rightarrow g_{x}(x,y)&=&g(x,y)-g(x-1,y)\\ \frac{\partial g(x,y)}{\partial y}\rightarrow g_{y}(x,y)&=&g(x,y)-g(x,y-1) \end{eqnarray*}

中心差分

前後のバランスを考えて偏微分を中心差分で近似する方法もある（【数値情報処理a】第6回 数値微分を参照）。 \begin{eqnarray*} \frac{\partial g(x,y)}{\partial x}\rightarrow g_{x}(x,y)&=&\frac{g(x+1,y)-g(x-1,y)}{2}\\ \frac{\partial g(x,y)}{\partial y}\rightarrow g_{y}(x,y)&=&\frac{g(x,y+1)-g(x,y-1)}{2} \end{eqnarray*} これを畳み込み演算で表すと、それぞれ以下のフィルタになる。 \[ F_{cx}= \left( \begin{array}{ccc} 0&0&0\\ -1/2&0&1/2\\ 0&0&0 \end{array} \right) = \frac{1}{2}\times \left( \begin{array}{ccc} 0&0&0\\ -1&0&1\\ 0&0&0 \end{array} \right) \] \[ F_{cy}= \left( \begin{array}{ccc} 0&-1/2&0\\ 0&0&0\\ 0&1/2&0 \end{array} \right) = \frac{1}{2} \times \left( \begin{array}{ccc} 0&-1&0\\ 0&0&0\\ 0&1&0 \end{array} \right) \]

分数（小数）があると実数の計算になり時間がかかるため、全体にかかる分数 $\frac{1}{2}$ を除いた整数の行列が実際の計算に使われる2)。

例

boats をグレースケール化して輪郭抽出を行う。

gradient_forward.R は一階偏微分を前進差分で近似して輪郭抽出を行うプログラムである。

gradient_forward.R
  2. # 輪郭抽出（前進差分）
  3.  
  4. library(imager)
  5.  
  6. # グレースケール画像
  7. im <- grayscale(boats)
  8.  
  9. # x方向の前進差分のフィルタ
  10. f.fx <- as.cimg(c(0, 0, 0, 0, -1, 1, 0, 0, 0), x = 3, y = 3)
  11. # y方向の前進差分のフィルタ
  12. f.fy <- as.cimg(c(0, 0, 0, 0, -1, 0, 0, 1, 0), x = 3, y = 3)
  13.  
  14. # 畳み込み演算
  15. im.fx = correlate(im, f.fx)
  16. im.fy = correlate(im, f.fy)
  17.  
  18. # 同時に4枚並べるレイアウト
  19. layout(t(1 : 4))
  20.  
  21. plot(im, main = "元の画像")
  22. plot(im.fx, main = "x方向の前進差分")
  23. plot(im.fy, main = "y方向の前進差分")
  24. plot(sqrt(im.fx ^ 2 + im.fy ^ 2), main = "輪郭の強さ")
  25.  
  26. # レイアウトを元に戻す
  27. layout(1)

$x$ 方向の前進差分と $y$ 方向の前進差分では輪郭がはっきりする方向が異なることが分かる。

小テスト②

Moodle Server（非公式）で第2回の小テスト②を受験しなさい。

Sobelフィルタ

一階微分による輪郭抽出では、画像に含まれるノイズも強調してしまう傾向にある。 ノイズを画像から取り除くには平滑化といわれる処理を行うとよいため、輪郭抽出とともに平滑化を同時に行うフィルタが考えられている。

そのようなフィルタの一つにSobel フィルタ がある。 畳み込み演算では以下のフィルタになる。 \[ F_{Sx}= \left( \begin{array}{ccc} -1/8&0&1/8\\ -1/4&0&1/4\\ -1/8&0&1/8 \end{array} \right) = \frac{1}{8}\times \left( \begin{array}{ccc} -1&0&1\\ -2&0&2\\ -1&0&1 \end{array} \right) \] \[ F_{Sy}= \left( \begin{array}{ccc} 1/8&1/4&1/8\\ 0&0&0\\ -1/8&-1/4&-1/8 \end{array} \right) = \frac{1}{8} \times \left( \begin{array}{ccc} 1&2&1\\ 0&0&0\\ -1&-2&-1 \end{array} \right) \]

全体にかかる分数 $\frac{1}{8}$ を除いた整数の行列を Sobel フィルタという。

二階微分による輪郭抽出

二度の微分を行うことを二階微分という。 二階微分は一階微分の変化の度合いを表す。 一階微分では山や谷になっているところを輪郭として抽出したが、そこは二階微分では正から負もしくは負から正へ $0$ を横切るところなので、そこを輪郭として抽出しても良い。

$g(x,y)$ の二階微分にも $x$ 方向と $y$ 方向があるので、やはり偏微分を行う。 $x,y$ の二階偏微分には四通りの組み合わせがある。 \begin{align*} \frac{\partial}{\partial x}\left(\frac{\partial g(x,y)}{\partial x}\right)&=\frac{\partial^{2}g(x,y)}{\partial x^{2}}\\ \frac{\partial}{\partial x}\left(\frac{\partial g(x,y)}{\partial y}\right)&=\frac{\partial^{2}g(x,y)}{\partial x\partial y}\\ \frac{\partial}{\partial y}\left(\frac{\partial g(x,y)}{\partial x}\right)&=\frac{\partial^{2}g(x,y)}{\partial y\partial x}\\ \frac{\partial}{\partial y}\left(\frac{\partial g(x,y)}{\partial y}\right)&=\frac{\partial^{2}g(x,y)}{\partial y^{2}} \end{align*}

これを行列で並べると \[ \left( \begin{array}{cc} \frac{\partial^{2}g(x,y)}{\partial x^{2}} & \frac{\partial^{2}g(x,y)}{\partial x\partial y} \\ \frac{\partial^{2}g(x,y)}{\partial y\partial x} & \frac{\partial^{2}g(x,y)}{\partial y^{2}} \end{array} \right) \] となる。 これをヘッセ行列という。 一般的な連続な関数では \[ \frac{\partial^{2}g(x,y)}{\partial x\partial y}=\frac{\partial^{2}g(x,y)}{\partial y\partial x} \] なので、ヘッセ行列の非対角成分は等しい。

デジタル画像では二階偏微分も差分で近似する。 例えば $\frac{\partial^{2}g(x,y)}{\partial x\partial y}$ は、$x$ 方向と $y$ 方向でそれぞれ前進差分を使うと以下になる。 \begin{eqnarray*} \frac{\partial^{2}g(x,y)}{\partial x\partial y}\rightarrow g_{xy}(x,y)&=&\{g(x+1,y+1)-g(x,y+1)\}\\ &&-\{g(x+1,y)-g(x,y)\}\\ &=&g(x,y)-g(x+1,y)\\ &&-g(x,y+1)+g(x+1,y+1) \end{eqnarray*} これを畳み込み演算で表すと、以下のフィルタなる。 \[ \left( \begin{array}{ccc} f(1,1)&f(2,1)&f(3,1)\\ f(1,2)&f(2,2)&f(3,2)\\ f(1,3)&f(2,3)&f(3,3) \end{array} \right) = \left( \begin{array}{ccc} 0&0&0\\ 0&1&-1\\ 0&-1&1 \end{array} \right) \]

ラプラシアン

ヘッセ行列は4通りあってそれぞれ方向が異なる。 そこで、ヘッセ行列の対角成分を足し合わせて使うことが多い。 \[ \frac{\partial^{2}g(x,y)}{\partial x^{2}}+\frac{\partial^{2}g(x,y)}{\partial y^{2}} =\left(\frac{\partial^{2}}{\partial x^{2}}+\frac{\partial^{2}}{\partial y^{2}}\right)g(x,y) =\Delta g(x,y) \] 演算子 $\Delta=\frac{\partial^{2}}{\partial x^{2}}+\frac{\partial^{2}}{\partial y^{2}}$ をラプラシアンという。

これを前進差分と後退差分を組み合わせて近似すると以下になる。 \begin{eqnarray*} \Delta\rightarrow L(x,y)&=&g(x,y-1)\\ &&+g(x-1,y)-4g(x,y)+g(x+1,y)\\ &&+g(x,y+1) \end{eqnarray*} 畳み込み演算で表すと、以下のフィルタになる。

\[ \left( \begin{array}{ccc} f(1,1)&f(2,1)&f(3,1)\\ f(1,2)&f(2,2)&f(3,2)\\ f(1,3)&f(2,3)&f(3,3) \end{array} \right) = \left( \begin{array}{ccc} 0&1&0\\ 1&-4&1\\ 0&1&0 \end{array} \right) \]

例

boats をグレースケール化して輪郭抽出を行う。

laplacian.R はラプラシアンで輪郭抽出を行うプログラムである。

laplacian.R
  2. # 輪郭抽出（ラプラシアン）
  3.  
  4. library(imager)
  5.  
  6. # グレースケール画像
  7. im <- grayscale(boats)
  8.  
  9. # ラプラシアンのフィルタ
  10. f.lap <- as.cimg(c(0, 1, 0, 1, -4, 1, 0, 1, 0), x = 3, y = 3)
  11.  
  12. # 畳み込み演算
  13. im.lap = correlate(im, f.lap)
  14.  
  15. # 同時に3枚並べるレイアウト
  16. layout(t(1 : 3))
  17.  
  18. plot(im, main = "元の画像")
  19. plot(im.lap, main = "ラプラシアン")
  20. plot(abs(im.lap), main = "ラプラシアン（絶対値）")
  21.  
  22. # レイアウトを元に戻す
  23. layout(1)

輪郭抽出の関数

imager には、あらかじめ輪郭抽出用の関数が用意されているので、実際はフィルタを作成しなくても輪郭抽出が簡単にできる。

関数 機能
imgradient(cimg変数) 一階偏微分の差分による畳み込み演算を行う。
オプション axes で偏微分の方向を指定できる（デフォルトでは axes = “xy”）。
　axes = “x” や axes = “y” だと $x$ 方向や $y$ 方向の結果を cimg 型で返す。
　axes = “xy” だと $x$ 方向と $y$ 方向の結果を同時に imlist 型で返す。
オプション scheme でグラディエントの方法を選択できる（デフォルトでは schme = 3）。
　1 は後退差分、2 は Sobel、3 は rotation invariant など。
imhessian(cimg変数) 二階偏微分の差分による畳み込み演算を行う。
オプション axes で偏微分の方向を指定できる（デフォルトでは axes = c(“xx”, “xy”, “yy”)）
　axes = c(“xx”, “xy”, “yy”) だと結果を imlist 型で返す。
imlap(cimg変数) ラプラシアンによる畳み込み演算を行う。

例

boats をグレースケール化して輪郭抽出を行う。

gradient.R は imgradient() で輪郭抽出を行うプログラムである。

gradient.R
  2. # 輪郭抽出（imgradient）
  3.  
  4. library(imager)
  5.  
  6. # グレースケール画像
  7. im <- grayscale(boats)
  8.  
  9. # グラディエント
  10. im.grad = imgradient(im)
  11.  
  12. plot(im.grad, layout = "row")
  13.  
  14. # plot(enorm(im.grad), main = "輪郭の強さ")
  15.  
  16. # もしくは
  17. # layout(t(1:2))
  18. # plot(im.grad$x, main = "x")
  19. # plot(im.grad$y, main = "y")
  20. # layout(1)

gradient_forward.R と同じことを imgradient() だけでできる。

課題

Moodle Server（非公式） で第2回の課題を行いなさい。

締め切り：2025年6月25日（水）20時

1)
カーネル、オペレータ、コンボリューション行列などとも呼ばれる。
2)
imager の場合は画素値が実数のため、あまり関係ない。