OpenCVとsklearnでk-means法を使う

# はじめに
前回の記事「[**OpenCVでトゥーンな加工をする**](https://qiita.com/mo256man/items/8a9bb7bb66e107fd7153)」のために減色について調べていたところ、k-means法による減色の記事を多数見つけることができた。k-meansといえば教師なし学習のクラスタリング手法ではないか。なぜ減色に使えるのだろう。
不思議に思った私はk-meansについて調べはじめた。

# k-means法について

## k-meansの考え方
紙と鉛筆で実感するキカガクの動画がわかりやすいので説明はそちらに任せることにする。

<iframe width="560" height="465" src="https://www.youtube.com/embed/8KhT_TCmLSw" frameborder="0" allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>

なるほど、k-meansは個々の点についてクラスタリングするだけでなく各クラスタの中心位置を得ることもできるわけだ。それが色における減色された代表色なのだな。はい、謎は解けた。終わり。

…というわけにはいかない。本番はここからだ。

## k-meansの特徴

 - クラスタ数はハイパーパラメータで人間様が決める必要がある。  
 自動でいい感じの値を算出する手法もいろいろ考案されてはいる。
 - 距離の二乗で計算するので、複雑な分布をクラスタリングすることはできない。
 - 距離の二乗和で評価するので、各ベクトルは標準化しておく必要がある（場合がある）。
 - 計算の第一手に依存するので部分最適に陥らないよう何度か計算する必要がある。


# Pythonでk-meansを使う
Pythonではscikit-learnやOpenCVが関数を持っている。
紙と鉛筆で作れるほどなので勉強のために関数をゼロから作っている人も少なくない。

## scikit-learnのk-means
scikit-learnではmodelを定義してfitするという機械学習でおなじみの使い方をする。
`sklearn.cluster.KMeans`はすべての引数にデフォ値が設定されているので省略しまくってお手軽に試すこともできる。クラスタ数が省略可能といっても自動で最適な値に決めてくれるわけではない。省略したら8になるだけだ。

```Python:
from sklearn.cluster import KMeans

model = KMeans(n_clusters=8, *, init='k-means++', n_init=10, 
               max_iter=300, tol=0.0001, verbose=0,
               random_state=None, copy_x=True, algorithm='auto')
```
`n_clusters=8`の次に`*`があるが、これは`*`以降はキーワード引数を指定する必要があるという意味。素直に2番目の引数として文字列ですらない記号`*`を指定するとエラーになる。というか、エラーになった。

　　参考　**[[python初心者向け]関数の引数のアスタリスク(*)の意味](https://dev.classmethod.jp/articles/what-does-asterisk-mean-at-args/)**

### 引数

  - n_clusters　クラスタ数。省略可能でデフォ値は8。
  - init　最初の一手の指定方法。省略可能でデフォ値は`'k-means++'`。
    - `'k-means++'`　kmeans++ の中心初期化手法。
    - `'random'`　　ランダムに選択される
    - numpy配列も指定できるが略
  - n_init　繰り返し回数。省略可能でデフォ値は10。
  - max_iter　繰り返しの最大回数。省略可能でデフォ値は300。
  - tol　toleranceの略で、相対許容誤差。省略可能でデフォ値は1e-4。
  - verbose　機械学習でよく見るオプション。デフォ値は0で、1にすると途中の詳しい計算を見ることができる。
  - random_state　デフォ値は0。ランダムシードを固定する場合は整数値を指定する。
  - copy_x　面倒なので[公式](https://scikit-learn.org/stable/modules/generated/sklearn.cluster.KMeans.html)を読んで。
  - algorithm　面倒なので[公式](https://scikit-learn.org/stable/modules/generated/sklearn.cluster.KMeans.html)を読んで。

### 使い方
`sklearn.cluster.KMeans()`でモデルを定義し、データ群を`fit()`させて使う。
`fit()`したモデルに`labels_`メソッドを使うことでラベルを取得するだけでなく、学習とクラスタリングを一度でおこなう`fit_predict()`もある。使いこなせば便利だとは思うが、普通に`fit()`を使うほうが良いだろう。

モデルだから訓練に使われていないデータを食わせてそのラベルを推定することもできる。これは後述の`cv2.kmeans()`にはない機能だ。
メソッドはほかにもあるが省略。

```Python:
model = KMeans(3)                      # クラスタ数=3とする
model.fit(X)                           # データXで訓練する

labels = model.labels_                 # 各データに付与されたラベル
centers = model.cluster_centers_       # 全n_clusters個のクラスタの中心位置
inertia = model.inertia_               # 各点とクラスタ中心の距離の二乗和
predict_labels = model.predict(new_X)  # 新たなデータのラベル
```

データ X および new_X は2次元配列。numpy配列でなくリストでも良いが、ベクトル数=1でも`[x1, x2, ...]`と単に列挙するのは不可で`[[x1], [x2], ...]`とする。
`model.labels_`および`model.predict()`は`(データ数,)`のshapeを持つ1次元のnumpy配列を返す。
`model.cluster_centers_`は`(クラスタ数, ベクトル数)`のshapeを持つ2次元のnumpy配列を返す。
`model.inertia_`は数値を返す。

## OpenCVのk-means

OpenCVでは以下のように使う。

```Python:
import cv2

compactness, labels, centers 
  = cv2.kmeans(data, K, bestLabels, criteria, attempts, flags[, centers] )
```

### 引数
 - data　データ。`np.float32`型のnumpy配列である必要がある。
 - K　　クラスタ数。
 - bestLabels　[チュートリアルで省略されている](http://labs.eecs.tottori-u.ac.jp/sd/Member/oyamada/OpenCV/html/py_tutorials/py_ml/py_kmeans/py_kmeans_opencv/py_kmeans_opencv.html)くらいなので割愛。使わないときはNoneとする。
 - criteria　繰り返し処理の終了条件。要素数3のタプル。精度のみ・回数のみの場合でも3個必要。
    - COUNT　繰り返し回数
      - cv2.TERM_CRITERIA_EPS　　　　　指定された精度に到達したら終了する
      - cv2.TERM_CRITERIA_MAX_ITER　　最大回数に到達したら終了する
      - cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER　　精度と回数の両方
   - MAX_ITER　繰り返しの最大回数
   - EPS　　　　精度
  - attempts　試行回数。
  - flags　　　最初の一手の指定方法。
    - cv2.KMEANS_RANDOM_CENTERS　ランダムに選択される
    - cv2.KMEANS_PP_CENTERS　　　kmeans++ の中心初期化手法が使われる
    - もう一つあるようだが略
  - centers　必須ではないので略。

### 出力
  - compactness　各点とクラスタ中心の距離の二乗和。  
  sklearnの`model.inertia_`に相当する。実際はこの値が使われることはほとんどないだろう。ではなぜこれが戻り値の先頭に位置しているかというと、これこそが繰り返し計算して最小にしたい値、すなわち損失関数の値だから。だと思う。
  - labels　　　　各点に付与された全K種類のラベル。  
 sklearnの`model.labels_`とは違い、`(データ数, 1)`のshapeを持つ2次元numpy配列となる。少々使いづらいが、元データが`(データ数, ベクトル数)`なのでそれと対応しているということで自分を納得させよう。
  - centers 　　　全K個のクラスタの中心位置。  
sklearnの`model.cluster_centers_`と同じく`(クラスタ数, ベクトル数)`のshapeとなる。


# コード
多くの解説サイトでは元データとしてアイリスのデータセットやskleanの`make_blobs`を使っているが、ここではそれらを使わず自分で作ることにする。こうすれば将来任意のデータを扱えるようになるだろう。

## サンプルデータを作るクラス
少し前なら関数としていたところだが、せっかく[クラスを習得した](https://qiita.com/mo256man/items/f81b7106184474e2a197)のでクラスで。

```Python:コード1
import numpy as np
from sklearn.cluster import KMeans
import random
import math
import cv2
import matplotlib.pyplot as plt

class SampleData():
    def __init__(self):
        self.n_clusters = 3
        self.width = 300
        self.height = 400
        self.colors_cv = [(0,0,255), (0,255,0), (255,0,0)]
        self.colors_skl = ["red", "green", "blue"]
        self.X = []

        # クラスタを3個作る
        self.random_plot(100, 100,  80, 100)
        self.random_plot(200, 100,  50, 100)
        self.random_plot(150, 300, 100, 200)

        # 単独の点をプロットする
        hazure_pts = [(25, 200), (275,200)]
        for pts in hazure_pts:
            self.X.append(pts)

        self.cnt = len(self.X)

    def random_plot(self, x0, y0, r0, cnt):
        for i in range(cnt):
            angle = 2 * math.pi * random.random()
            r = r0 * random.random()
            x = x0 + r * math.cos(angle)
            y = y0 + r * math.sin(angle)
            self.X.append((x, y))
```

これで作られるサンプルデータは**リスト**だ。

## `sklearn.cluster.KMeans()`でクラスタリング

リストXを直接`KMeans()`に食わせている。
`matplotlib.pyplot`でグラフ化するにあたりxの値のリストとyの値のリストをあらためて作っている。データ群Xの個々の要素をxとしたかったが、xはほかで使うのでelementの略でelmとした。こういうときどういう変数名にすればいいんでしょ。

```Python:コード2
def sklearn_kmeans(data):
    model = KMeans(data.n_clusters)
    model.fit(data.X)
    labels = model.labels_
    centers = model.cluster_centers_
    colors = [data.colors_skl[i] for i in labels]

    # 以下、グラフ化
    fig, ax = plt.subplots()
    ax.set_xlim(0, data.width)
    ax.set_ylim(0, data.height)
    ax.set_aspect("equal")
    ax.invert_yaxis()                   # OpenCV画像と合わせるためy軸を反転する

    # 各点をプロットする
    x = [elm[0] for elm in data.X]
    y = [elm[1] for elm in data.X]
    ax.scatter(x, y, s=5, c=colors)

    # 各クラスタの中心点をプロットする
    cx = [elm[0] for elm in centers]
    cy = [elm[1] for elm in centers]
    ax.scatter(cx, cy, s=100, c="black", marker="o")

    # 各クラスタの中心点の座標を表記する
    for i in range(data.n_clusters):
        msg = f"({cx[i]:.2f},{cy[i]:.2f})"
        ax.annotate(msg, (cx[i], cy[i]))

    # モデルを元に未知の点をクラスタリングする
    new_X = [(100, 200), (200, 200)]
    predict_labels = model.predict(new_X)
    nx = [elm[0] for elm in new_X]
    ny = [elm[1] for elm in new_X]
    ncolors = [data.colors_skl[i] for i in predict_labels]
    ax.scatter(nx, ny, s=100, c=ncolors, marker="D")

    plt.show()
```

## `cv2.kmeans()`でクラスタリング

こちらではデータを`np.float32`型のnumpy配列に直している。
関数の引数はあまりこだわっていないのであしからず。
`cv2.kmeans()`で得られたクラスタリング結果を`matplotlib.pyplot`で表現することはもちろん可能だが、せっかくなのでここではOpenCV画像にしている。

```Python:コード3
def cv2_kmeans(data):
    X = np.array(data.X, np.float32)
    criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 10, 1.0)
    _, labels, centers = cv2.kmeans(X, data.n_clusters, bestLabels=None, criteria=criteria,
                                    attempts=10, flags=cv2.KMEANS_RANDOM_CENTERS)

    # 以下、画像化
    h, w = data.height, data.width
    image = np.full((h, w, 3), (255,255,255), np.uint8)

    # 各点をプロットする
    for elm, label in zip(X, labels):
        x, y = int(elm[0]), int(elm[1])
        color = data.colors_cv[label[0]]       # labelではなくlabel[0]
        cv2.circle(image, (x,y), 2, color, -1)

    # 各クラスタの中心点をプロットする
    for i in range(data.n_clusters):
        cx, cy = int(centers[i][0]), int(centers[i][1])
        color = data.colors_cv[i]
        cv2.circle(image, (cx,cy), 5, (0,0,0), -1)
        msg = f"({centers[i][0]:.2f},{centers[i][1]:.2f})"
        cv2.putText(image, msg, (cx,cy), cv2.FONT_HERSHEY_DUPLEX, 0.5, (0,0,0))

    cv2.imshow("cv2_kmeans", image)
    cv2.waitKey(0)
    cv2.destroyAllWindows()
```

## 全体コード

クラス化・関数化しているのでそれらを合体するのみ。

```Python:
import numpy as np
from sklearn.cluster import KMeans
import random
import math
import cv2
import matplotlib.pyplot as plt
# 以上、コード1で示したやつ

class SampleData():
    # コード1で示したやつ

def sklearn_kmeans(data):
    # コード2で示したやつ

def cv2_kmeans(data):
    # コード3で示したやつ

if __name__ == "__main__":
    data = SampleData()
    sklearn_kmeans(data)  # sklearnによるクラスタリング
    cv2_kmeans(data)      # OpenCVによるクラスタリング
```

## 実行結果
赤、緑、青がラベル0,1,2を意味するが、最初の一手がランダムなので実行するたびにラベルが変わる。区分けはするがどれがどのラベルに属するかは決められない。これが教師なし学習だ。

外れ値として(x,y)=(25, 200)および(275,200)の点を置いたが、k-means法ではこれらも計算によっていずれかのラベルが与えられている。今回は詳しく述べなないが、[DBSCAN](https://ja.wikipedia.org/wiki/DBSCAN)ならばハイパーパラメータの指定により外れ値であると結論づけることができる。

sklearnのみ、学習には使われていない新たな点(x,y)=(100, 200)および(200, 200)もクラスタリングされている。直感的に納得できる結果かどうかは別として。

|sklearn|OpenCV|
|:-:|:-:|
|![sklearn_kmeans.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/520216/0958dcc5-9b23-226d-a12b-1a5f3076b0d0.png)|![cv2_kmeans.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/520216/192a2262-9f4c-3b0f-cc1a-d2786fb55421.png)|

# 終わりに
クラスタリングしてmatplotlib.pyplotのグラフとOpenCV画像を作るにあたり、二次元データをそのまま使ったりデータ型を変換したりx・yの各ベクトルに分けて使ったり(x1,y1),(x2,y2)と個々のデータを使ったりとさまざまな加工をおこなった。
この経験が重要だ。この苦労をしないと、ネットに転がっているサンプルコードを動かすことはできても自分が分析したいデータに応用することができないという羽目に陥ってしまう。

…とかっこいいことを書いたところで今回はここまで。減色に応用する前に~~力尽きてしまったので~~けっこうなボリュームになったので。