線形判別分析を用いた画像の2値化

17456 ワード
statistics Python 判別分析 画像処理 2値化 Python テキストリンク

本記事のゴール

フィッシャーの線形判別分析 (Fisher's linear discriminant analysis) を理解し、それを用いた画像の2値化を行えるようになる。

使用する環境

ソフトウェア名 バージョン
Python 3.4 or 3.5
Pillow 3.1
Numpy 1.10
Matplotlib 1.5.0

フィッシャーの線形判別分析

概要

  • データをクラス分類する線形判別式のパラメタを求める手法。
  • 線形判別式は $y=w^{T}x$ で表される。ここで $y$ はクラス、$w$ はパラメタ、$x$ はデータ。
  • $w$ は「クラス間平均の距離を大きくし、クラス内分散を小さくする」ように求められる。
  • ベイズ最適分類器とも深い関係がある。
  • 画像の2値化の他、データの次元削減などにも応用できる。

目的関数

下記関数を最大化するようなパラメタ $w$ を見つける。

\begin{equation}
J(w) = \frac{w^{T}S_{B}w}{w^{T}S_{W}w}
\end{equation}
  • $w$: 線形判別式のパラメタ。
  • $S_{W}$: クラス内共分散行列。
  • $S_{B}$: クラス間共分散行列。
  • $T$: 転置。

クラス内共分散行列

$N$ 個のデータを各クラスのデータ数が $N_{i} (i=0,...,k-1)$ となるように$k$個のクラスに分類した時の、各クラスの共分散行列の平均。

\begin{equation}
S_{W}=\sum_{i=0}^{k-1}\frac{N_{i}}{N}S_{i}
\end{equation}
  • $S_{i}$: クラス $i$ の共分散行列。

クラス間共分散行列

各クラスの平均ベクトルの共分散行列。

\begin{equation}
S_{B}=\sum_{i=0}^{k-1}\frac{N_{i}}{N}(\mu_{i}-\mu)^{T}(\mu_{i}-\mu)
\end{equation}
  • $\mu_{i}$: クラス $i$ の平均ベクトル。
  • $\mu$: データ全体の平均ベクトル。

データ全ての共分散行列、クラス内共分散行列、クラス間共分散行列の関係

クラス内共分散行列 ($S_{W}$) とクラス間共分散行列 ($S_{W}$) の和は、データ全ての共分散行列 ($S$) と等しい。以降、データ全ての共分散行列を全共分散行列と呼ぶ。

\begin{align}
S_{W} + S_{B} &=\sum_{i=0}^{k-1}\frac{N_{i}}{N}S_{i}+\sum_{i=0}^{k-1}\frac{N_{i}}{N}(\mu_{i}-\mu)^{T}(\mu_{i}-\mu)\\
&=\frac{1}{N}\sum_{i=0}^{k-1}(N_{i}(\sum_{j \in C_{i}}\frac{x_{j}^{T}x_{j}}{N_{i}}-\mu_{i}^{T}\mu_{i})+N_{i}\mu_{i}^{T}\mu_{i})-\mu^{T}\mu\\
&=\frac{1}{N}\sum_{i=0}^{k-1}(\sum_{j \in C_{i}}x_{j}^{T}x_{j})-\mu^{T}\mu\\
&=\frac{1}{N}\sum_{l=0}^{N-1}x_{l}^{T}x_{l}-\mu^{T}\mu\\
&=S
\end{align}
  • $C_{i}$: i番目のクラス (に含まれるデータの集合)。
  • $x_{j}$: j番目のデータ。

ラグランジアンによるパラメタの決定

目的関数 $J(w)=\frac{w^{T}S_{B}w}{w^{T}S_{W}w}$ の $w$ を $\alpha w$ に置き換えてもその値は変化しないので、$w^{T}S_{W}w=1$ と置き換え可能であり、$J(w)$ を最大化するパラメタ $w$ を求める問題は、$w^{T}S_{W}w=1$ という制約のもとで $w^{T}S_{B}w$ を最大化する $w$ を求める問題に帰着する。

\begin{equation}
argmax_{w}(w^{T}S_{B}w),
\hspace{3mm}subject\hspace{3mm}to\hspace{3mm}
w^{T}S_{W}w=1
\end{equation}

これを最小化問題に変換する。

\begin{equation}
argmin_{w}(-\frac{1}{2}w^{T}S_{B}w),
\hspace{3mm}subject\hspace{3mm}to\hspace{3mm}
w^{T}S_{W}w=1
\end{equation}

そうすると、ラグランジアンは次式のようになる。

\begin{equation}
L=-\frac{1}{2}w^{T}S_{B}w+\frac{1}{2}\lambda (w^{T}S_{W}w - 1)
\end{equation}

$w$ で $L$ を微分し、それを $0$ とおけば次式が得られる。これを一般化固有値問題と呼ぶ。

\begin{equation}
S_{B}w=\lambda S_{W}w
\end{equation}

$S_{W}$ を正則と仮定すると、両辺に左からその逆行列 $S_{W}^{-1}$ をかけて通常の固有値問題を得る。

\begin{equation}
S_{W}^{-1}S_{B}w=\lambda w
\end{equation}

一般化固有値問題を解いた結果得られる固有値のうちの最大のものに対応する固有ベクトルを $w$ として採用する。この選択基準は、本最適化問題の双対問題を考えることによって得られる。

フィッシャーの判別分析を用いた画像の2値化

概要

  • グレースケール画像を白黒画像に変換。
  • 閾値 $t$ を設定し、$t$ より画素値が小さいものを黒、大きいものを白に分類。
  • 閾値 $t$ をフィッシャーの線形判別分析に基づいて決定する。そうすると、その意味で最適な閾値を自動的に選択可能。
  • この方法は、元画像が本来白黒で表現されるような画像を想定。

目的関数

前章のフィッシャーの判別式の目的関数 $J(w)$ は、グレースケール画像のデータが1次元であることに注意しながら、クラス内分散、クラス間分散、全分散の関係を使うと次式を得る。ここで、パラメタ $w$ は1次元であることに注意する。

\begin{equation}
J(w) = \frac{\sigma_{B}}{\sigma_{W}}=\frac{\sigma_{B}}{\sigma -\sigma_{B}}
\end{equation}
  • $\sigma_{W}$: クラス内分散。
  • $\sigma_{B}$: クラス間分散。
  • $\sigma$: 全分散。

閾値の決定

全分散 $\sigma$ は一定であるので、$J(w)$ を最大化するような $\sigma_{B}$ を構成する閾値 $t$ を決定すれば良い。

$\sigma_{B}$ は次式で求めることができる。

\begin{equation}
\sigma_{B}=\frac{N_{0}}{N}(\mu_{0}-\mu)^2 + \frac{N_{1}}{N}(\mu_{1}-\mu)^{2}
\end{equation}
  • $N$: 画像全体の画素数
  • $N_{0}$: 黒として分類される画素の数
  • $N_{1}$: 白として分類される画素の数
  • $\mu$: 画像全体の画素値の平均
  • $\mu_{0}$: 黒として分類される画素の画素値の平均
  • $\mu_{1}$: 白として分類される画素の画素値の平均

グレースケール画像の画素値の最小と最大をそれぞれ $x_{min}$、$x_{max}$ とすれば、$x_{min} \leq t \leq x_{max}$ である。このような $t$ 全てに対して $\sigma_{B}$ を計算し、$J(w)$ を最大化するものを見つけ閾値とすれば良い。

プログラム

# This file is part of "Junya's self learning project about digital image processing."
#
# "Junya's self learning project about digital image processing"
# is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# "Junya's self learning project about digital image processing"
# is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with Foobar.  If not, see <http://www.gnu.org/licenses/>.
#
# (c) Junya Kaneko <[email protected]>

from PIL import Image
import numpy as np
from matplotlib import pyplot as plt


def discriminant_analysis(image):
    _image = image.reshape(image.shape[0] * image.shape[1])
    ave = np.average(_image)
    t = j0 = 0
    for _t in range(np.min(_image) + 1, np.max(_image) + 1):
        blacks = _image[_image < _t]
        whites = _image[_image >= _t]
        sigma_b = len(blacks) * np.power(np.average(blacks) - ave, 2) /len(_image) + len(whites) * np.power((np.average(whites) - ave), 2) / len(_image)
        j1 = sigma_b / (np.var(_image) - sigma_b)
        if j0 < j1:
            t = _t
            j0 = j1
    return t


if __name__ == '__main__':
    image = np.array(Image.open('img/test.jpg').convert('L'), dtype=np.uint8)
    bin_image = np.array(image)
    t = discriminant_analysis(image)

    for i in range(image.shape[0]):
        for j in range(image.shape[1]):
            bin_image[i, j] = 255 if image[i, j] >= t else 0

    plt.subplot(1, 2, 1)
    plt.title('Original')
    plt.imshow(image, cmap='Greys_r')

    plt.subplot(1, 2, 2)
    plt.title('Binary')
    plt.imshow(bin_image, cmap='Greys_r')

    plt.show()

ソースコードは Github にある。

実行結果は下図のようになる。元画像 (Original) は 字典 Wiki のものを使用した。

参考文献