import cv2
import numpy as np
import matplotlib.pyplot as plt
from OpenCV.plot_img import cv2plot
img = cv2.imread("Lenna.png")
fig, ax = plt.subplots(figsize=(6,6))
ax = cv2plot(img, ax=ax)
ax.set_title("Sample Image: Lenna")
plt.show()
2. ディジタル画像の撮影
2-4-1. グレースケール画像
- 元々グレースケールの画像
- カラー画像データはより多くの情報を含むが、カラー情報が必要ない事もある。(生物分野だと必要ない事が多いらしい。)
- 単一チャネルでのアルゴリズムが既に確立している場合
- 計算時間を短縮したい場合
- 特定チャンネルのデータに興味がある場合(蛍光など)
- RGBの3チャネルの情報をまとめてグレースケール変換する手法
- 単純:(R+G+B)/2
- 人間の光受容を考慮:(2R+4G+B)/7
opencv_ch = ["red", "green", "blue"]
fig, axes = plt.subplots(ncols=3, nrows=1, figsize=(18,6))
for i,(ax,ch) in enumerate(zip(axes,opencv_ch)):
ax.imshow(img[:,:,i], cmap='gray')
ax.set_title(f"{ch} channel.", fontsize=20, color=ch)
plt.tight_layout()
plt.show()
cv2_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
pure_gray = np.mean(img, axis=2, dtype=np.int)
cie_gray = np.average(img, axis=2, weights=[0.298912,0.586611,0.114478]).astype(np.int)
grays = [cv2_gray, pure_gray, cie_gray]
titles = ["OpenCV default", "Average", "Weighted Average"]
fig, axes = plt.subplots(ncols=3, nrows=1, figsize=(18,6))
for i,ax in enumerate(axes):
ax.imshow(grays[i], cmap='gray')
ax.set_title(titles[i], fontsize=20)
plt.tight_layout()
plt.show()
2-4-2. 標本化と量子化
- 標本化(sampling):アナログ信号から、離散的な位置におけるアナログ値を取り出す処理。
- 量子化(quantization):標本値を有限分解能の数値に変換する処理。(ex.8ビット量子化→(0-255),1ビット量子化→2値画像)
- 標本化により画像は一般的に横 $M$ 個・縦 $N$ 個の画素(pixel)で表現される。
- 原点 $(0,0)$ は"左上"である事が多い。
- 3次元ではさらに奥行き方向(voxel)にも標本化を行う。
- 動画ではさらに時間軸方向(temporal)にも標本化を行う。
- 画像の標本化では、2次元的な像(アナログ画像)の光強度に対して、縦横に等間隔の格子状に配置した標本点(sampling point)における光強度を電圧値として取り出す。
- 標本点の間隔を標本化間隔(sampling interval)と呼び、標本点におけるアナログ値を標本値(sampling value)と呼ぶ。
- 標本化間隔が周波数成分の半分よりも大きいと、エリアシングノイズ(aliasing noise)が発生する事がある。
3. 画像の性質と色空間
3-1-1. 画像の統計量
- 横軸に画素値、縦軸にそれぞれの画素値の頻度を取る濃淡ヒストグラム(単にヒストグラム)が一般的。
- その他、平均値・分散・中央値(median)・最頻値(mode)など。
- ヒストグラムが詰まっている → コントラストが低い
fig, (axgray, ax3ch) = plt.subplots(ncols=2, nrows=1, figsize=(12,4), sharey=True, sharex=True)
axgray.hist(np.ravel(cv2_gray), bins=256, color="black", density=True)
for i,ch in enumerate(opencv_ch):
ax3ch.hist(np.ravel(img[:,:,i]), bins=256, color=ch, density=True, alpha=0.3)
axgray.set_ylabel("Frequency", fontsize=16)
axgray.set_title("Grayscale Image."), ax3ch.set_title("RGB 3ch Image.")
plt.show()
3-1-2. コントラストとシャープネス
- コントラストは、画像のヒストグラムの分布の広がりを表す。画素値の最大値を $I_{\max}$、最小値を $I_{\min}$ とするとき、以下のような求め方がある。
$$
\begin{aligned}
C_1 &= I_{\max} - I_{\min}\\
C_2 &= \left(I_{\max} - I_{\min}\right) / \left(I_{\max} + I_{\min}\right)\\
C_3 &= I_{\max} / I_{\min}
\end{aligned}
$$
- どの求め方も外れ値に対して敏感である。そのため、頻度が一定値よりも小さな画素値を考慮しない、などの工夫が必要になる。
- $C_2$ を使うと、$I_{\min}=0$ のとき、$\forall I_{\max}C_2=1$ となるので注意が必要である。
- シャープネス(鮮鋭度;sharpness)は、画像を見た時に感じる鮮鋭感を表す尺度である
- 鮮鋭度が高い画像:エッジ付近の画素値変化が急激で、画像中の細やかな部分まで鮮明に観察する事ができる。
- 鮮鋭度が低い画像:ピントが合ってないような印象を受け、エッジ付近の画素値変化が緩やかで、画像中の細やかな部分を読み取りにくい。
4. 画素ごとの濃淡変換
4-1. 明るさ・コントラストの変換
ディジタル画像の各画素は、その濃淡を表す値(画素値)を持っている。そこで、画像の濃淡を変化させるためには、「入力画像のそれぞれの画素値」に対し、「出力画像の画素値」をどのように対応づけるか指定すれば良い。
そのような対応関係を与える関数のことを階調変換関数(gray-level transformation function)、また、それをグラフで表したものをトーンカーブ(tone curve)と呼ぶ。
pixels = np.arange(0,2**8,1)
colors = ["black", "green", "red", "blue"]
gammas = [1,1.5,2,3]
gamma_transform = lambda x,gamma: 255*(x/255)**gamma
plt.figure(figsize=(6,6))
for gamma,color in zip(gammas, colors):
plt.plot(pixels, gamma_transform(pixels, 1/gamma), color=color, linestyle=":", label=f"$\gamma={1/gamma:.2f}$")
plt.plot(pixels, gamma_transform(pixels, gamma), color=color, label=f"$\gamma={gamma:.1f}$")
plt.legend(), plt.title("Gamma tone Curve")
plt.show()
有名なトーンカーブに、上図のガンマ補正(ガンマ変換)がある。これは、以下の形で表される。
$$y = 255\left(\frac{x}{255}\right)^{\gamma}$$ここで、$y=x$ よりも下側($\gamma > 1$)だと出力画像が入力画像に比べて明るくなり、上側($\gamma < 1$)はその逆が起こる事がわかる。
また、コントラストが低い画像はS字トーンカーブを用いて変換される事がある。これは、binの密度が高い部分を引き延ばす事で、コントラストを上げている。(ヒストグラム平坦化(histogram equalization)と似ている。)
4-3-2. 疑似カラー
- R,G,B各チャンネルに対して適当に異なるトーンカーブを用いることで、グレースケール画像に対して疑似的な色(疑似カラー(pseudo color))をつける事ができる。
- グレースケール画像の解釈を容易にする。(区別のつきにくい微妙な濃淡の違いを視覚的に表現できる。)
- 正解画像データの作成時に利用する。
- 診断補助に用いる。
def apply_pseudo_color(x):
r = 255 if x>=192 else 0 if x<=128 else 4*(x-128)
g = 255 if 64<=x<=192 else 4*x if x<=64 else 4*(255-x)
b = 255 if x<=64 else 0 if x>=128 else 4*(128-x)
return (r,g,b)
R,G,B = np.vectorize(apply_pseudo_color)(pixels)
plt.plot(R, color="red")
plt.plot(G, color="green")
plt.plot(B, color="blue")
plt.xlabel("Input pixel", fontsize=16), plt.ylabel("Output pixel", fontsize=16), plt.title("Simple Pseudo Color Transformer.", fontsize=14)
plt.xticks([0,255]),plt.yticks([0,255])
plt.show()
pseudo_img = np.stack(np.vectorize(apply_pseudo_color)(cv2_gray), axis=-1)
fig, axes = plt.subplots(ncols=2, nrows=1, figsize=(12,6))
axes[0].imshow(cv2_gray, cmap='gray'), axes[0].set_title("Grayscale", fontsize=16)
axes[1].imshow(pseudo_img), axes[1].set_title("Pseudo Color", fontsize=16)
plt.tight_layout()
plt.show()
4-4-1. 画像間演算
複数の入力画像($f_1,f_2$)の同じ画素に対して AND
,OR
などの論理演算(logic operation)や四則演算(arithmetic operation)などの演算を用いて出力画像($g$)を生成する。例としては、
- アルファブレンディング(alpha blending): $$g = \alpha f_1 + \left(1-\alpha\right)f_2$$
- ディゾルブ(dissolve)・オーバーラップ(overlap):$\alpha$ の値を画素の位置で変化させる。
- エンボス(emboss):濃淡を反転した画像($\bar{f_1}$)を平行移動させ($\bar{f_1}_{x+dx}$)元の画像と足し合わせる。 $$g = f_1 + \bar{f_1}_{x+dx} - 128$$
def emboss(f1, dx=0.01, constant_values=0):
"""
@params f1: shape=(H,W) must be the grayscale image.
@params dx: (int,float) distance of translation.
"""
H,W = f1.shape
if dx<1:
dx = min(H,W)*dx
dx = int(dx)
f1_inv = 255-f1
f2 = np.pad(f1_inv, (0,dx), 'constant', constant_values=constant_values)[-H:,-W:]
g = np.clip(f1+f2-128, 0, 255).astype(int)
return g
embossed_img = emboss(cv2_gray, dx=2)
fig, axes = plt.subplots(ncols=2, nrows=1, figsize=(12,6))
axes[0].imshow(cv2_gray, cmap='gray'), axes[0].set_title("Grayscale", fontsize=16)
axes[1].imshow(embossed_img, cmap='gray'), axes[1].set_title("Embossed", fontsize=16)
plt.tight_layout()
plt.show()
5. 領域に基づく濃淡変換(空間変換)
5-1. 空間フィルタリング
- 「画素ごと(pixel-to-pixel operation)に」画素値を変換するのではなく、「周囲の画素の画素値も利用して(area-to-pixel operation)」画素値を変換する。
- 大別すると、「線形フィルタ(linear filter)」と「非線型フィルタ(nonlinear filter)」がある。
- 線形フィルタは以下で表される。 $$g(i,j) = \sum_{n=-W}^W\sum_{m=-W}^Wf(i+m,j+n)h(m,n)$$
- 非線型フィルタは、上記以外で表されるもの全てである。
def apply_filter(img, f):
"""
@params img: (ndarray) Grayscale Image, shape=(W,H)
@params f : (ndarray) Filter: shape=(FW,FH)
@params out: (ndarray) Output Image, shape=(OW,OH)
"""
W,H = img.shape
FW,FH = f.shape
OW,OH = (W-FW+1,H-FH+1)
out = np.zeros(shape=(OW,OH))
for i in range(OW):
for j in range(OH):
out[i,j] = np.sum(f*img[i:i+FW,j:j+FH])
out = np.clip(out, 0, 255).astype(int)
return out
5-2. 平滑化
平滑化(smoothing)を行うと、画像の濃淡変化を滑らかにする事ができる。したがって、画像に含まれるノイズなどの不要な濃淡変化を軽減するために用いられる。
def mk_average_filter(filter_size):
"""filter_size: (int)"""
if filter_size<=0:
return np.ones(shape=(1,1))
return np.full(shape=(filter_size, filter_size), fill_value=1/filter_size**2)
filter_sizes = [0,3,5]
fig, axes = plt.subplots(ncols=3, nrows=1, figsize=(18,6))
for ax,fs in zip(axes, filter_sizes):
f = mk_average_filter(fs)
ax.imshow(apply_filter(cv2_gray, f), cmap='gray')
ax.set_title(f"Filter Size: {fs}", fontsize=20)
plt.tight_layout()
plt.show()
※ だんだんとボケているのがわかる。
他にも、
- 加重平均フィルタ(weighted average filter)
- ガウシアンフィルタ(gaussian filter)
- メディアンフィルタ(median filter)
- 特定方向の平滑化
などがある。
5-3. エッジ抽出
エッジ抽出(edge extraction)は、画像中で明るさが急激に変化するエッジ部分を取り出す事であり、
- 微分フィルタ:連続関数の場合は $$f^{\prime}(x) = \lim_{h\rightarrow0}\frac{f(x+h)-f(x)}{h}$$ で微分が計算できるが、デジタル画像の場合には、注目画素とその隣接画素との差分で置き換えられる。
- プリューウィットフィルタ(Prewitt filter): 微分フィルタでは、エッジ部分を抽出すると同時に「画像に含まれるノイズを強調する」傾向がある。そこで、ノイズを抑えながらエッジを強調するために、「 ある方向への微分」と「それに直交する方向に関する平滑化」を施すことを考えたフィルタ。
- ソーベルフィルタ(Sobel filter): プリューウィットフィルタ(Prewitt filter)における「縦方向への平滑化」のときに、中央に重みを付けた平均化を行う方法。
- 2次微分フィルタ(second derivative filter): 上記の微分フィルタの概念を拡張し、微分を2回繰り返したフィルタ。
- ラプラシアンフィルタ(Laplacian filter): 2次微分の値を用いてラプラシアンを $$\frac{\partial^2}{\partial x^2}f(x,y) + \frac{\partial^2}{\partial y^2}f(x,y)$$ 求める事ができるので、これを用いて方向に依存しないエッジを直接得る。
- LoGフィルタ: 一般に、ラプラシアンは微分を繰り返すため、ノイズを強調してしまう。そこで、まずガウシアンフィルタを適用してある程度の平滑化を行ったのち、ラプラシアンフィルタを施す事がよく行われる。2次元ガウス分布のラプラシアンは $$h_{log}(x,y) = \frac{x^2+y^2-2\sigma^2}{2\pi\sigma^6}\exp\left(-\frac{x^2+y^2}{2\sigma^2}\right)$$ で得られるので、この $h_{log}(x,y)$ を係数とするフィルタを用意すれば良い。
prewitt_filter_h = 1/6*np.array([[-1,0,1],[-1,0,1],[-1,0,1]])
prewitt_filter_v = np.rot90(prewitt_filter_h)
sobel_filter_h = 1/8*np.array([[-1,0,1],[-2,0,2],[-1,0,1]])
sobel_filter_v = np.rot90(sobel_filter_h)
filters = [
[prewitt_filter_h,prewitt_filter_v],
[sobel_filter_h ,sobel_filter_v]
]
titles = [
["Prewitt filter (Vertical)", "Prewitt filter (Horizontal)"],
["Sobel filterr (Vertical)", "Sobel filter (Horizontal)"],
]
fig, axes = plt.subplots(ncols=2, nrows=2, figsize=(8,8))
for i,ax_cols in enumerate(axes):
for j,ax in enumerate(ax_cols):
ax.imshow(apply_filter(cv2_gray, filters[i][j]), cmap='gray')
ax.set_title(f"Filter: {titles[i][j]}")
plt.tight_layout()
plt.show()
derivative_filter = np.array([[0,0,0],[-1/2,0,1/2],[0,0,0]])
second_derivative_filter = np.array([[0,0,0],[1, -2, 1],[0,0,0]])
laplacian_filter = np.array([[0,1,0],[1, -4, 1],[0,1,0]])
filters = [derivative_filter, second_derivative_filter, laplacian_filter]
titles = ["Derivative Filter (Horizontal)", "Second Derivative Filter (Horizontal)", "Laplacian Filter"]
fig, axes = plt.subplots(ncols=3, nrows=1, figsize=(18,6))
for i,ax in enumerate(axes):
ax.imshow(apply_filter(cv2_gray, filters[i]), cmap='gray')
ax.set_title(f"Filter: {titles[i]}", fontsize=20)
plt.tight_layout()
plt.show()
5-4. 鮮鋭化
上のフィルタは、画像のエッジ部分だけを取り出すようなものであった。ここでは、元の画像の濃淡を残したままエッジを強調する、画像の鮮鋭化(sharpening)処理を考える。
- 入力画像に対して平滑化処理を施し、その結果を元の画像から引く → 元の画像のエッジ部分が取り出されたような画像が得られる。
- それを元の入力画像と足し合わせることにより、画像のエッジが強調された(鮮鋭化された)画像が得られる。
このような処理をアンシャープマスキング(unsharp masking)と呼ぶ。
なお、この二段階の処理は上記のエッジ抽出のフィルター同様一つのフィルタ(鮮鋭化フィルタ(sharping filter))で表現する事ができる。
def mk_sharpening_filter(k, filter_size=3):
f = np.full(shape=(3,3), fill_value=-k/9)
f[1][1] = 1+8/9*k
return f
ks = [4,9,18]
fig, axes = plt.subplots(ncols=3, nrows=1, figsize=(18,6))
for ax,k in zip(axes, ks):
f = mk_sharpening_filter(k)
ax.imshow(apply_filter(cv2_gray, f), cmap='gray')
ax.set_title(f"Filter Size: {fs}", fontsize=20)
plt.tight_layout()
plt.show()
6. 周波数領域におけるフィルタリング
画像のフーリエ変換
2次元画像は、2次元空間で定義される関数 $f(x,y)$ として表されるので、2次元フーリエ変換を行う事で、
- $x,y$ で表される空間領域(spatial domain)
- $u,v$ で表される周波数領域(frequency domain)
を関連づける事ができる。
$$ \begin{cases} \begin{aligned} F(u,v) &= \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} f(x,y)\exp\left\{-j2\pi\left(ux+vy\right)\right\}dxdy & \left(\text{Fourie transform}\right)\\ f(x,y) &= \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} F(u,v)\exp\left\{-j2\pi\left(ux+vy\right)\right\}dudv & \left(\text{inverse Fourie transform}\right)\\ \end{aligned} \end{cases}\\ j = \sqrt{-1} $$なお、上は連続関数のフーリエ変換であるが、ディ離散的フーリエ変換(discrete Fourier transform)ジタル画像を変換する場合には、離散的逆フーリエ変換(inverse discrete Fourier transform)を用いる。
さらに、その計算を高速に実行する高速フーリエ変換(Fast Fourier transform; FFT)があり、実用的によく用いられている。
なお、お気持ち的には、
「2次元の様々な周波数の三角関数(ex.$z = \sin(ax)\cdot\sin(by)$)」が入力にどの程度含まれているかを周波数ごとに分解する事で調べる。
周波数フィルタリング
- 画像はフーリエ変換により、周波数領域における別の形で表現する事ができる。
- フーリエ変換の結果は、画像に含まれるそれぞれの周波数成分を示している。
- フーリエ変換の各周波数成分の大きさを各成分ごとに変える事により、元の画像の性質を変化させる事ができる。
このような処理を周波数フィルタリング(frequency filtering)と呼ぶ。
- 元の画像をフーリエ変換したものを $F(u,v)$
- フィルタリングの出力を $G(u,v)$
- 周波数フィルタを $H(u,v)$
とするとき、周波数フィルタリングは
$$G(u,v) = F(u,v)H(u,v)$$として表される。なお、周波数フィルタリングでは、入力の各周波数成分 $F(u,v)$ が、周波数ごとにフィルタ $H(u,v)$ と掛け算され、出力 $G(u,v)$ となる。
周波数フィルタの具体例として、「画像に含まれる空間周波数成分のうち、ある特定の範囲の周波数成分を残す」ようなタイプのフィルタとして、
- ローパスフィルタ(lowpass filter):低周波数成分のみ残す。
- ハイパスフィルタ(highpass filter):高周波数成分のみ残す。
- バンドパスフィルタ(bandpass filter):ある中間的な周波数の範囲を残す。
- ガウス分布型のローパスフィルタ:空間領域のガウシアンフィルタに対応。
などがある。一方、上記のハイパスフィルタでは、画像の直流成分を含む低周波数成分を除去してしまうため、画像の平均的な明るさが保たれないという問題がある。
そこで、画像の低周波数成分はそのまま保ちつつ、高周波数成分を強調するフィルタとして高域強調フィルタ(high-emphasis filter)がある。
9. 2値画像処理
9-1. 2値化
目的に応じて具体的な手法は様々。
判別分析法(discriminant analysis method)(大津の2値化)が有名。
変数 | 意味 |
---|---|
$$t$$ | しきい値 |
$$m_1,\sigma_1^2,\omega_1$$ | 黒画素クラスの平均・分散・画素数 |
$$m_2,\sigma_2^2,\omega_2$$ | 白画素クラスの平均・分散・画素数 |
$$\sigma_{\omega}^2 = \frac{\omega_1\sigma_1^2 + \omega_2\sigma_2^2}{\omega_1+\omega_2}$$ | クラス内分散(within-class variance) |
$$\begin{aligned}\sigma_b^2 &= \frac{\omega_1(m_1-m_t)^2 + \omega_2(m_2-m_t)^2}{\omega_1 + \omega_2}\\&= \frac{\omega_1\omega_2(m_1-m_2)^2}{(\omega_1+\omega_2)^2}\end{aligned}$$ | クラス間分散(between-class variance) |
$$\sigma_2^2=\sigma_b^2 + \sigma_{\omega}^2$$ | 全分散(total variance) |
$$\frac{\sigma_b^2}{\sigma_w^2} = \frac{\sigma_b^2}{\sigma_t^2 - \sigma_b^2}$$ | 分離度(separation metrics) |
この定義の下で、分離度を最大にするしきい値 $t$ を求めたいが、このとき、しきい値 $t$ に関係なく全分散 $\sigma_t^2$ は一定のため、クラス間分散 $\sigma_b^2$ が最大になるような $t$ を決めれば良い。
def maximize_separation_metrics(img):
thresholds = np.arange(0,256,1)
between_class_variances = np.zeros_like(thresholds)
N = len(img)
for i,t in enumerate(thresholds):
black_img = img[img>=t]
white_img = img[img<t]
Nb=len(black_img); Nw=len(white_img)
if Nb==0 or Nw==0: continue
Mb=np.mean(black_img); Mw=np.mean(white_img);
between_class_variances[i] = Nb*Nw*(Mb-Mw)**2/(N)**2
return thresholds,between_class_variances
thre,bcv = maximize_separation_metrics(cv2_gray)
opt_i = np.argmax(bcv)
opt_t = thre[opt_i]
opt_bcv = bcv[opt_i]
plt.plot(thre, bcv, color="black")
plt.scatter(opt_t, opt_bcv, color="red", s=100, label="Optimal")
plt.title("Discriminant Analysis Method"),plt.xlabel("Thresholds"), plt.ylabel("between-class variance")
plt.show()
binary_img = np.where(cv2_gray>=opt_t,255,0)
fig, axes = plt.subplots(ncols=2, nrows=1, figsize=(12,6))
axes[0].imshow(cv2_gray, cmap="gray"), axes[0].set_title("Grayscale", fontsize=16)
axes[1].imshow(binary_img, cmap="gray"), axes[1].set_title("Optimal Binary Image.", fontsize=16)
plt.show()
9-2-1. 連結性
- 2値画像では、値1を持つひとまとまりの領域を定義するため、連結(connection)という概念を用いる。
- ディジタル画像では、注目画素に対して上下左右の画素を4近傍(4-connected neighbor)と呼ぶ。
- 4近傍に斜め方向の近傍を加えたものを8近傍(8-connected neighbor)と呼び、その近傍に対して注目の画素の連結を定義したものを8連結(8-connection)と呼ぶ。
- 連結している画素の集合を連結成分(connected component)と呼ぶ。
- 対象の連結成分の中にあり、背景に連結していない白画素の集合を穴(hole)と呼ぶ。
9-2-2. 輪郭追跡
- 連結成分の境界を求めることを輪郭追跡(contour tracking)と呼ぶ。
- 詳しいアルゴリズムに関しては
- 3S-Bioinformatics-Programming/c++/README.mdに説明がある。
- kadai16.cppにC++のコードがある。
9-2-4. 収縮・膨張処理
- 背景または穴に接する対象の画素をひとまわりはぎとる処理を、収縮(erosion)と呼ぶ。この処理は、連結性を保存しない。
- 背景または穴に接する対象の画素に、画素をひとまわり加える処理を、膨張(dilation)と呼ぶ。この処理も連結性を保存しない。
- 同じ回数だけ、膨張して収縮する処理をクロージング(closing)と呼ぶ。これを施すことで画像の小さな穴を除くことができる。
- 同じ回数だけ、収縮したのち膨張する処理をオープニング(opening)と呼ぶ。これを施すことで画像の小さな連結成分を除くことができる。
9-2-6. ラベリング
- ラスタキャノンによるラベリング
- 輪郭追跡による距離変換画像(distance transform image)
9-3. 線画像のベクトル化
- 連結成分の連結性を保存したまま画素を削る処理を繰り返すことを細線化(thinning)と呼ぶ。
- しきい値処理をして得られた2値画像の線分は、一般に幅を持っており、線分の特徴を抽出しにくい。
- 細線化を行い、線幅$1$の2値画像に変換する。細線化された2値画像の画素は、特徴点である以下の3種類に分けられる。
- 端点(end point)
- 分岐点(branch point)
- 通過点(passing point)
- 幾何学的特徴を見るために、「端点-端点」・「端点-分岐点」・「分岐点-分岐点」というように、細線を画素列に分割し、さらにその画素列を直線近似し、始点と終点の座標値で表現するベクトル型の数値データに変換する。
10. 領域処理
10-2-3. ミーンシフトを用いた領域分割処理
- ミーンシフト(mean shift)法は、画像処理の分割や対象画像の追跡に用いられる手法である。
- カーネル密度推定法に求められる密度分布の勾配を計算することで重心の移動方向を決定する。
- アルゴリズムは以下:
- カラー画像の$N$個の各画素の位置$\mathbf{x}_i$とその画素値を$\boldsymbol{\nu}_i=(H_i,S_i,I_i)$とし、画素位置-画素値を結合した5次元空間内の点 $\mathbf{z}_i=(\mathbf{x}_i,\boldsymbol{\nu_i})$ を考えてみて、ミーンシフト法で各画素をクラスタリングする。
- 画素位置に関するバンド幅$h_s$、画素値に関するバンド幅$h_r$を与え、全ての$\mathbf{z}_i$にミーンシフトを行い、収束位置$\mathbf{z}^c=(\mathbf{x}^c,\boldsymbol{\nu}^c)$を計算する。
- $\mathbf{x}_i$ の画素値を、収束位置の値($H^c,S^c,I^c$)に置き換えることによって画像の平滑化ができる。
- なお、ミーンシフトは、$\mathbf{x}^s,\mathbf{x}^r$ をそれぞれ5次元ベクトル(位置を表す2次元+色を表す3次元)の空間に対応するものとし、カーネル密度推定を以下の通り定義する。($c$ は全体の積分値を$1$にするための係数である。) $$f(\mathbf{x})= \frac{c}{Nh_s^2h_r^3}\sum_{i=1}^Nk\left(\left|\frac{\mathbf{x}^s-\mathbf{x}^s_i}{h_s}\right|^2\right)k\left(\left|\frac{\mathbf{x}^r-\mathbf{x}^r_i}{h_r}\right|^2\right)$$
10-2-4. 対象物と背景の間のエッジを利用した領域分割処理
- 対象領域がエッジで囲まれているときは、エッジを閉曲線として抽出するスネーク(snakes)と呼ばれる手法によって領域分割が行える。
- 画像中の対象領域の多くは、背景との境界において画素値の差を生じ、それがエッジとなって現れる。
- 必然的に対象物体の領域は、エッジを境界線とした閉領域になる。
- 対象物体を囲む閉曲線を初期値として与え、徐々に閉曲線が縮んでいき、エッジの境界に沿って貼り付くことによって、対象物体の領域を抽出する。
11. パターン・図形・特徴の検出とマッチング
11-3-1. コーナー検出
ハリスのコーナー検出(Harris corner detector)は、画像からコーナーを検出する代表的な手法である。
- 入力画像$I$に対して、$x$軸方向の勾配画像$I_{x}$、$y$軸方向の勾配画像$I_{y}$を生成する。ここでは、ガウス関数$G(\sigma)$を$x,y$の各方向で微分した$G_x(\sigma),G_y(\sigma)$を画像に畳み込むことで、勾配画像を求める。 $$I_x = G_x(\sigma)\ast I,\quad I_y =G_y(\sigma)\ast I$$
- 各勾配画像の積により、各方向における勾配の大きさを算出する。 $$I_{x2}=I_x\cdot I_x,\quad I_{2y}=I_y\cdot I_y,\quad I_{xy}=I_x\cdot I_y$$
- $I_{x2},I_{y2},I_{xy}$ の局所領域における勾配の総和$S_{x2},S_{xy},S_{y2}$を求める。$S_{x2},S_{y2},S_{xy}$ は局所領域における単純な総和でも良いが、以下に示すようなガウス関数$G(\sigma^{\prime})$による重み付き和を用いることが多い $$S_{x2} = G(\sigma^{\prime})\ast I_{x2},\quad S_{y2} = G(\sigma^{\prime})\ast I_{y2},\quad S_{xy} = G(\sigma^{\prime})\ast I_{xy}$$
- 画素$(x,y)$における局所領域の勾配の総和$S_{x2}(x,y),S_{y2}(x,y),S_{xy}(x,y)$を要素にもつ行列 $\mathbf{M}(x,y)$ を定義する。 $$\mathbf{M}(x,y)=\left[\begin{array}{rrr} S_{x2}(x,y) & S_{xy}(x,y) \\ S_{xy}(x,y) & S_{y2}(x,y) \\ \end{array}\right]$$
- 画素$(x,y)$がコーナーである場合、行列$\mathbf{M}$の固有値$\lambda_1,\lambda_2$は共に大きな値となる。そこで、コーナー関数 $R$を以下のように定義する。なお、$k$は調整パラメータであり、$0.04\sim0.06$ が最適値とされている。 $$\begin{aligned} R&=\operatorname{det}\mathbf{M}-k\left(\operatorname{tr}\mathbf{M}\right)^2\\ \operatorname{det}\mathbf{M} &= \lambda_1\cdot\lambda_2,\quad\operatorname{tr}\mathbf{M} = \lambda_1+\lambda_2 \end{aligned}$$
- 上式の値が局所的な最大値となる画素をコーナーとして検出する。実際には、画像中に非常に多くの局所的な最大値が存在するため、適当なしきい値を設け、有効なコーナーだけを選択する。また、類似度補間手法のように$2$次関数を当てはめることで、コーナーのサブピクセル位置を推定することもできる。
11-3-2. DoG画像を用いた特徴点とスケールの検出
画像中に拡大縮小があると、画像間の特徴点領域の濃淡パターンが変化するため、特徴点の対応づけができない。
そこで、特徴点とその領域の大きさを表すスケールを検出する必要があり、複数のDoG(Difference-of-Gaussian)画像を用いて計算することができる。
DoGはLoG(Laplacian-of-Gaussian)を近似したものであり、スケールの異なるガウス関数$G(\sigma)$と入力画像を畳み込んだ平滑化画像 $L$ の差分により、DoG画像 $D(\sigma)$ を求める。
$$ \begin{aligned} D(\sigma) &= \left(G(k\sigma)-G(\sigma)\right)\ast I\\&= L(k\sigma)-L(\sigma)\\ L(\sigma) &= G(\sigma)\ast I \end{aligned} $$ここで、$k$は$\sigma$の増加率であり、スケールを少しずつ大きくして複数のDoG画像を求める。
11-4-1. スケールと回転に不変な特徴記述(SIFT)
SIFT特徴記述では、回転に不変な特徴量を記述するためにオリエンテーションを求める。
特徴量記述の際にオリエンテーションにより向きの正規化をすることで、回転に不変な特徴量を得ることができる。
13. 動画像処理
13-1. 差分画像を用いた移動物体検出
- 動画像の情報抽出処理では、異なる時間に撮影された$2$枚の画像の差を観察することによって、変化情報を得ることができる。
- 2毎の画像において同じ位置にある画素値の差の絶対値を画像としたものを差分画像(subtraction image)と呼ぶ。
- 背景差分法(background subtraction method)
- 移動物体がない状態の画像を固定カメラで背景画像として取り込む。
- 移動物体が入った画像から1の画像を差分し、移動物体の領域が0以外の値を持った差分画像を得る。
- 得られた差分画像の画素値に対してしきい値処理を行い、2値画像を得る。
- クロージングおよびオープニングを用いてこれらを取り除き、移動物体の領域を得る。
- 得られた2値画像を利用して対象物体の領域に位置する画素を取り出すことにより、移動物体の画像を得る。
- 移動物体がない理想的な背景画像を得ることができない場合、移動物体を撮影した異なる時間の3枚の画像を用いて、フレーム間差分法(frame subtraction method)によって移動物体領域を取り出すことができる。
- AとB、BとCの差分画像を作成し、しきい値処理を施し2値画像ABとBCを得る。
- 2値画像ABとBCの論理積処理(AND)を行い、ABとBCの共通領域を取り出す。
- 観測される対象シーンの中に樹木などの定常的に変動している物体が写り込んでしまう場合がある。そのような場合、今までの手法だと樹木の小枝も移動物体として検出されてしまう。そこで、画素値の定常的な変動を考慮して移動物体を検出する統計的背景差分法(statistical background subtraction method)がある。
13-2. オプティカルフロー
- 異なる時間に撮影された2枚の画像を用いて、静止カメラで撮影した移動物体の運動を解析することができる。
- 異なる時間に撮影された2枚の画像間での対象の移動量をベクトルデータとして表現したものをオプティカルフロー(optical flow)と呼ぶ。
- ブロックマッチング法(block matching method)は、テンプレートマッチングを用いてオプティカルフローを求める方法である。
- 勾配法(gradient-based method)は、連続する2枚の画像での対象物の移動量が微小であることを前提にオプティカルフローを求める手法である。
16. 画像符
16-3-1. ハフマン符号化
ハフマン符号化(Huffman coding)は、以下の手順によって処理される。
- 出現確率の最も小さい2つのシンボルを選択する。
- 出現確率の大きい方に符号0、小さい方に符号1を割り当て、部分木を作成する。
- 2つのシンボルにおける出現確率の和を出現確率とする新たなシンボルに結合する。
- 1~3を繰り返す。
16-5-1. ランレングス符号化
2値画像において、画素値を横方向にスキャンすると、0の連続と1の連続が繰り返すことがわかる。
これを、同じ画素値が連続したもの(ラン)と、その連続数(レングス)の組の並びで符号化することをランレングス符号化(run-length coding)と呼ぶ。