PythonでZhang-Suenの細線化アルゴリズムを実装

**注意**：結構この記事読んでくれる方が多いんですが、こちらの記事に書いてあるコードだと激烈に遅いので、アルゴリズムの説明を読み終わったら[こっち](https://qiita.com/hrs1985/items/c97eaeb1bc9e4a75dd39)のコード使ってください。

細線化は、(主には2値化した)画像を幅1ピクセルの線画像に変換する操作です。
用途は色々あるのですが、画像処理ライブラリのOpenCVになぜか細線化のメソッドが無かったので、
勉強も兼ねて自分で実装しました。
速さが欲しい場合[こちら](http://qiita.com/hrs1985/items/c97eaeb1bc9e4a75dd39)へどうぞ。

#細線化のアルゴリズム
細線化にはいくつかのアルゴリズムが提案されています。

1. 田村のアルゴリズム
2. Zhang-Suenのアルゴリズム
3. Nagendraprasad-Wang-Guptaのアルゴリズム

それぞれ癖があるのですが、
今回は2番のZhang-Suenのアルゴリズムを使いました。
(ロジックが簡単で書きやすそうだったので)

#Zhang-Suenのアルゴリズム
細線化に用いる画像は、あらかじめ2値化(黒を1、白を0)にしておきます。
##前準備
###注目している(P1)ピクセルの周囲のピクセルに番号を振る
[P9][P2][P3]
[P8][P1][P4]
[P7][P6][P5]

###2つの関数を定義する
f1(P1): P2,P3,P4,P5,P6,P7,P8,P9,P2.と並べて順番に見ていったとき、0の次が1となっている場所の個数
f2(P1): P2~P9の中の1の個数

##処理
###ステップ1
以下の5つの条件をすべて満たすピクセルを記録しておく。

1. 黒である
2. f1(P1)がちょうど1
3. f2(P1)が2以上6以下
4. P2, P4, P6のいずれかが白
5. P4, P6, P8のいずれかが白

記録したピクセルを白に変える。
###ステップ2
以下の5つの条件をすべて満たすピクセルを記録しておく。

1. 黒である
2. f1(P1)がちょうど1
3. f2(P1)が2以上6以下
4. P2, P4, __P8__のいずれかが白
5. __P2__, P6, P8のいずれかが白

記録したピクセルを白に変える。

以上の2つのステップを、
どちらのステップでも変更点が無くなるまでくり返す。

#Pythonでの実装
使用する関数を作っておきます。

```functions.py
# Zhang-Suenのアルゴリズムを用いて2値化画像を細線化します
def Zhang_Suen_thinning(binary_image):
    # オリジナルの画像をコピー
    image_thinned = binary_image.copy()
    # 初期化します。この値は次のwhile文の中で除かれます。
    changing_1 = changing_2 = [1]
    while changing_1 or changing_2:
        # ステップ1
        changing_1 = []
        rows, columns = image_thinned.shape
        for x in range(1, rows - 1):
            for y in range(1, columns -1):
                p2, p3, p4, p5, p6, p7, p8, p9 = neighbour_points = neighbours(x, y, image_thinned)
                if (image_thinned[x][y] == 1 and
                    2 <= sum(neighbour_points) <= 6 and # 条件2
                    count_transition(neighbour_points) == 1 and # 条件3
                    p2 * p4 * p6 == 0 and # 条件4
                    p4 * p6 * p8 == 0): # 条件5
                    changing_1.append((x,y))
        for x, y in changing_1:
            image_thinned[x][y] = 0
        # ステップ2
        changing_2 = []
        for x in range(1, rows - 1):
            for y in range(1, columns -1):
                p2, p3, p4, p5, p6, p7, p8, p9 = neighbour_points = neighbours(x, y, image_thinned)
                if (image_thinned[x][y] == 1 and
                    2 <= sum(neighbour_points) <= 6 and # 条件2
                    count_transition(neighbour_points) == 1 and # 条件3
                    p2 * p4 * p8 == 0 and # 条件4
                    p2 * p6 * p8 == 0): # 条件5
                    changing_2.append((x,y))
        for x, y in changing_2:
            image_thinned[x][y] = 0        
    
    return image_thinned

# 2値画像の黒を1、白を0とするように変換するメソッドです
def black_one(binary):
    bool_image = binary.astype(bool)
    inv_bool_image = ~bool_image
    return inv_bool_image.astype(int)

# 画像の外周を0で埋めるメソッドです
def padding_zeros(image):
    import numpy as np
    m,n = np.shape(image)
    padded_image = np.zeros((m+2,n+2))
    padded_image[1:-1,1:-1] = image
    return padded_image

# 外周1行1列を除くメソッドです。
def unpadding(image):
    return image[1:-1, 1:-1]

# 指定されたピクセルの周囲のピクセルを取得するメソッドです
def neighbours(x, y, image):
    return [image[x-1][y], image[x-1][y+1], image[x][y+1], image[x+1][y+1], # 2, 3, 4, 5
             image[x+1][y], image[x+1][y-1], image[x][y-1], image[x-1][y-1]] # 6, 7, 8, 9

# 0→1の変化の回数を数えるメソッドです
def count_transition(neighbours):
    neighbours += neighbours[:1]
    return sum( (n1, n2) == (0, 1) for n1, n2 in zip(neighbours, neighbours[1:]) )

# 黒を1、白を0とする画像を、2値画像に戻すメソッドです
def inv_black_one(inv_bool_image):
    bool_image = ~inv_bool_image.astype(bool)
    return bool_image.astype(int) * 255
```

コードの本体はここから

```main.py
import matplotlib.pyplot as plt
import cv2
import numpy as np

# 画像を読み込みます
image = cv2.imread('image.jpg')
# グレースケールに変換します
image_gray = cv2.cvtColor(image, cv2.COLOR_RGB2GRAY)
# ガウシアンフィルタをかけます
blur = cv2.GaussianBlur(image_gray,(5,5), 3)
# 大津のアルゴリズムで2値化します
ret,th2 = cv2.threshold(blur,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)
# 2値化画像の黒を1、白を0に変換します。外周を0で埋めておきます。
th2 = padding_zeros(th2)
new_image = black_one(th2)
# Zhang Suenアルゴリズムによる細線化を行います
result_image = inv_black_one(Zhang_Suen_thinning(new_image))
new_image = inv_black_one(unpadding(new_image))
# 結果を出力します
plt.subplot(121), plt.imshow(image_gray, cmap='gray')
plt.title("input image")
plt.subplot(122), plt.imshow(result_image, cmap='gray')
plt.title("thinned image")
plt.show()
```

#結果
処理結果です。
NやQに少し怪しげな部分が残ってますが、大体できているようです。
![figure_1.png](https://qiita-image-store.s3.amazonaws.com/0/134222/c2769ad8-0725-90f2-60e1-c1d0e4812f07.png)

最初は自力で実装したのですが、
実装後に見つけた[githubに上がっていたコード](https://github.com/linbojin/Skeletonization-by-Zhang-Suen-Thinning-Algorithm "Zhang-Suen Thinning Algorithm, Python and Matlab Implementation")の方が綺麗だったので倣って修正を加えました。（消滅していました。2016年11月26日確認。リンクを修正しました。消滅していませんでした。2023年2月16日）

2018年7月17日追記
mainの中で細線化前の画像のpaddingと細線化後のunpaddingが抜けていたので追加しました。コメントにてご指摘いただいたDo_you_1istenさんありがとうございます。