[C/C++]Zhang-Suenアルゴリズムによる細線化の実装[コード付]

どーも
今回の画像処理は前回までと毛色を変えて細線化を実装してみたいと思います。
細線化とは？からコードの実装までやっていくのでぜひ最初から最後までご覧ください！

前提

前提として僕はXcodeでコードを書いて、実行しています。
opencvはhomebrewからダウンロードして、xcodeからパスをつないだり必要なファイルなどをブチ込んだり、、、とバンバンに使える状態から始めています。
Xcodeにopencvを導入する方法はgoogkeで調べればいくらでもあるので、まだ導入できてない人はまずそこをお願いします。（気が向いたらこのサイトでもopencvの導入の記事を書くかも

細線化とは？

まずはここから解説しましょう。
細線化とはズバリその名前の通り線を細くすることです。
具体的には、、

こんな感じ。白背景に描かれたアルファベットの「A」の線が細くなっていますよね。
これをC++で実装していきたいと思います。
ちなみにこの細線化がどういう時に役立つかというと、

ラベリング
面積算出

などに役立ちます。
というのもあるイラストに対してラベリングや面積を計算する際に、輪郭線が太いとその輪郭線自体も領域を持つことになるので上の操作の精密度が下がりますよね。そんな時に細線化の処理をかけることで輪郭線を１pxにすることで解決できます。

Zhang-Suenアルゴリズム

今回の記事ではZhang-Suenアルゴリズムを用いて細線化を行っていきます。
ここでは何でこのような操作で線が細くなるのかは解説しません（気になる人は自分で調べて）
事実としてつらつら書いていきます。。

Zhang-Suenアルゴリズムの具体的な手順

画素の取り扱い

ある注目画素に対して、以下の画像のように近傍画素に番号を振っていきます

関数A(P1)

画素P1に対して、
P2,P3,,,P9,P2と順番に画素をみた時に、値が0から1に変化した時(白の画素の次が黒だった)回数をカウントしてその値を返す

関数B(P1)

画素P1に対して、
P2からP9のうち画素の値が1（つまり黒色）のものの数をカウントしてその値を返す

処理α

細線化を行う画像に対してラスタ走査（x=0,y=0から順に走査）を行っていき、注目している画素が次の条件を全て満たしていればその位置を記録する。

P1=1（注目画素が黒）
A(P1)=1
2<=B(P1)<=6
P2×P4×P6=0
P4×P6×P8=0

すべての画素に対してチェックを行った後、記録された位置の画素を全て白に変更する。

処理β

P1=1（注目画素が黒）
A(P1)=1
2<=B(P1)<=6
P2×P4×P8=0
P2×P6×P8=0

すべての画素に対してチェックを行った後、記録された位置の画素を全て白に変更する。

メインの処理

上記の処理αと処理βを交互に行なっていく。
処理α→処理βの連続で記録された画素がない場合（つまり連続で画素の変更がなかった場合）全体の処理を終了する。（この時細線化が完了している）

コードの実装

画像の読み書き

// カラー画像を読み込む
void Read_Color( const char *name, unsigned char img[HEIGHT][WIDTH][3] )
{
    cv::Mat readImage( HEIGHT, WIDTH, CV_8UC3);
    readImage = cv::imread(name, 1);

    for( int w=0; w<WIDTH; w++ )
        for( int h=0; h<HEIGHT; h++ )
            for( int c=0; c<3; c++ )
                img[h][w][c] = readImage.data[(w+h*WIDTH)*3+c];
}

// カラー画像をファイルへ保存
void Write_Color( const char *fname, unsigned char img[HEIGHT][WIDTH][3] )
{
    cv::Mat IMAGE( HEIGHT, WIDTH, CV_8UC3 );

    for( int w=0; w<WIDTH; w++ )
        for( int h=0; h<HEIGHT; h++ )
            for( int c=0; c<3; c++ )
                IMAGE.data[(w+h*WIDTH)*3+c] = img[h][w][c];
    
    cv::imwrite( fname , IMAGE);
}

void showImage(unsigned char img[HEIGHT][WIDTH]){
    cv::Mat showImage( HEIGHT, WIDTH, CV_8U );

    for( int w=0; w<WIDTH; w++ ){
        for( int h=0; h<HEIGHT; h++ )
        {
            showImage.data[w+h*WIDTH] = img[h][w];
        }
    }
    
    cv::imshow("aaa", showImage);
    cv::waitKey(0);
}

こちらの関数で不明な点があればこちらの記事をどうぞ！↓

二値化

void binary(unsigned char org[HEIGHT][WIDTH], unsigned char byn[HEIGHT][WIDTH]){
    for(int w=0; w<WIDTH; w++){
        for(int h=0; h<HEIGHT; h++){
            if((int)org[h][w] < 200)
                byn[h][w] = 0;
            else
                byn[h][w] = 255;
        }
    }
}

ここらへんは前回までの画像処理の記事で書いたので解説は割愛。

アルゴリズムの説明の時点でそうでしたが、
今回は二値化画像（白黒）で考えていきます。

近傍8画素をまとめる関数

まず後の処理をスムーズに行うために近傍8画素を一つの配列の中に入れていきましょう。

上の画像の矢印の順番に入れていくと関数A()が記述しやすくなるので

void f_cre(int y, int x, unsigned char img[HEIGHT][WIDTH], int f[9]){

    f[0] = img[y - 1][x];
    f[1] = img[y - 1][x + 1];
    f[2] = img[y][x + 1];
    f[3] = img[y + 1][x + 1];
    f[4] = img[y + 1][x];
    f[5] = img[y + 1][x - 1];
    f[6] = img[y][x - 1];
    f[7] = img[y - 1][x - 1];
    f[8] = img[y - 1][x];
}

このようにしました。
のちに説明しますが、画像でのP2の画素の値をf[0]とf[8]の2つに格納しているのがミソです。

関数A()

上で紹介した関数ですね。とりあえずコード

int A_P(int f[9]){
    int count = 0;
    for (int i = 0; i < 8; i++) {
        if ((f[i+1]-f[i]) == -255)
            count++;
    }
    return count;
}

引数のint f[9]は上の関数で作成した近傍８画素を格納する配列です。
関数Aの処理の手順に従って白から黒に変わる回数を数えていくのですが、fの中には矢印の順番に画素値が格納されているので

for (int i = 0; i < 8; i++) {
        if ((f[i + 1] - f[i]) == -255)
            count++;
    }

このように簡単にfor文で書けちゃいます。
中身は一つ先の画素値との差を計算して-255ならばカウントを増やす感じです。これは黒(0)引く白(255)の値が-255だからですね。（他の３パターンでは-255にならない）
またここでf[8]にf[0]の値を入れておくことでP2とP9の比較も上のfor文で解決できます。

関数B()

とりあえずコード

int B_P(int f[9]){
    int count = 0;
    for (int i = 0; i < 8; i++) {
        if (f[i] == 0)
            count++;
    }
    //cout<<"B:"<<count<<endl;
    return count;
}

こっちはもっと簡単ですね。
同じようにf[9]を受け取って黒(0)の画素を数えて返しているだけです。

記録した位置の画素を変更する関数

こちらは処理αと処理βの最後に行う画素変更の際に用いる関数です。

void change(unsigned char img[HEIGHT][WIDTH], bool reco[HEIGHT][WIDTH]) {
    for (int w = 1; w < (WIDTH - 1); w++)
        for (int h = 1; h < (HEIGHT - 1); h++)
            if (reco[h][w])
                img[h][w] = 0;
}

位置記録用配列recoの中身に沿って、画像格納配列imgの値を変えていきます。

処理αの条件判定

処理αの関数を書く前に条件を判定する関数を作成しましょう。
5つの条件のうち、2〜５番を判定する関数にします。（１はあえてここでは判定しない）

A(P1)=1
2<=B(P1)<=6
P2×P4×P6=0
P4×P6×P8=0

enum p_num{
    p2,p3,p4,p5,p6,p7,p8,p9//記述をわかりやすくするための列挙型
};

bool judge_alpha(int y,int x,unsigned char img[HEIGHT][WIDTH]){
    int f[9]={0};
    f_cre(y, x, img, f);
    int A = A_P(f);
    int B = B_P(f);

    if((A==1)&&(2<=B)&&(B<=6)&&((f[p2]*f[p4]*f[p6])==0)&&((f[p4]*f[p6]*f[p8])==0))
        return true;
    
    return false;
}

まずは上のenumについて。
処理αと処理βの条件4,5はどちらも近傍画素を取り出して比較するものです。ただ画像の矢印のように配列fに値を格納したので、配列の添字と画像での添字が一致していないのでこのまま記述すると後から手直ししにくです。なので列挙型でp2〜p9までを0〜７に対応するように定義してこの関数内では用いています。
列挙型についてわからない方はこちらの記事をどうぞ！

では関数について。
引数で受け取っているy,xは注目画素の座標でimgは画像が格納されている配列です。
とりあえず上で説明した関数を用いて注目がそのA()とB()を計算します。
それらの値と条件の計算式を用いてif文の条件を設定。全て満たしていればtrue、そうでなければfalseが返るようなbool型の関数にします。
ここで先ほどのenumで定義した列挙子を用いて条件を書いています。

処理α

とりあえずコードです。

int alpha(unsigned char img[HEIGHT][WIDTH]) {
    bool reco[HEIGHT][WIDTH] = { false };
    int flag = 0;

    for (int h = 1; h < (HEIGHT - 1); h++){
        for (int w = 1; w < (WIDTH - 1); w++) {
            if((int)img[h][w]==0){
                if(judge_alpha(h, w, img)){
                    flag++;//画素の変更が起きた時にインクリメント
                    reco[h][w] = true;//画素の変更が起きる場所を記録
                }
            }
        }
    }
    change(img, reco);
    return flag;
}

引数で受け取っているimgは画像の配列です。
まず変更する画素の位置を記録する用のbool型配列recoと処理全体で画素の変更が起きたかを判定するための変数flagを定義します。

中の二重ループがメインの処理ですね。
wとhはそれぞれ画像の縦と横を示していてラスタ走査をこのループで行なっているわけです。
ここでwとhの範囲なのですが

for (int w = 1; w < (WIDTH - 1); w++)
        for (int h = 1; h < (HEIGHT - 1); h++)

このように画像に対して一回り小さい範囲のみのラスタ操作となっています。
なぜこのようにしたかというと画像の辺や角の画素は近傍８画素が定義できないからです。
なので一回り小さい範囲のみを走査しています。

次はループの中身です

if(img[h][w]){
    if(judge_alpha(h, w, img)){
          flag++;//画素の変更が起きた時にインクリメント
          reco[h][w] = true;//画素の変更が起きる場所を記録
     }
}

まず最初のifで注目画素が黒か白かを判定しています。
黒なら1なのでif文の条件を満たすということですね。ここの処理が説明で定義した条件１にあたります。
なんで処理αの条件判定の関数の中にこの処理を入れなかったかというと、計算量を削減するためです。注目画素が黒か白かは近傍8画素とは関係ありませんし、関数A()、B()の値とも関係ないですよね。なので簡単な条件１だけ外に出して先に判定させることで無駄な計算を省略しているわけです。

そして条件１を通ったら中のif文で処理αの条件判定の関数を評価して全ての条件を満たしているか判定します。
もし満たしているならばflagをインクリメントし、記録用の配列recoをtrueにして処理を終了します。

そして最後です

change(img, reco);
return flag;

記録した位置の画素を変更していきます。その際には関数changeにて処理します。
また返り値はflagとします。
何も変更がなければ0、そうでなければ1以上の値が返るわけですね。

処理βの条件判定

ここら辺はαと同様なので説明は割愛します

bool judge_beta(int y,int x,unsigned char img[HEIGHT][WIDTH]){
    int f[9]={0};
    f_cre(y, x, img, f);
    int A = A_P(f);
    int B = B_P(f);

    if((A==1)&&(2<=B)&&(B<=6)&&((f[p2]*f[p4]*f[p8])==0)&&((f[p2]*f[p6]*f[p8])==0))
        return true;
    
    return false;
}

処理β

int beta(unsigned char img[HEIGHT][WIDTH]) {
    bool reco[HEIGHT][WIDTH] = { false };
    int flag = 0;

    for (int h = 1; h < (HEIGHT - 1); h++){
        for (int w = 1; w < (WIDTH - 1); w++) {
            if((int)img[h][w]==0){
                if(judge_beta(h, w, img)){
                    flag++;//画素の変更が起きた時にインクリメント
                    reco[h][w] = true;//画素の変更が起きる場所を記録
                }
            }
        }
    }
    change(img, reco);
    return flag;
}

main関数

int main(int argc, const char * argv[]) {
    unsigned char img[HEIGHT][WIDTH]={0};
    cv::Mat saveImage( HEIGHT, WIDTH, CV_8U );
    readImage("thinning_5.png", img);
    
    binary(img, img);
    writeImage("gray_thin.png", img);

    
    int flag = 1;
    while (flag) {
        flag = 0;
        int a = alpha(img);//変更があったら1以上になる
        int b = beta(img);//変更があったら1以上になる
        flag = a + b;//どちらにも変更がない場合のみ０
        cout<<flag<<endl;
        showImage(img);
    }
    
    writeImage("thinning.png", img);
    return 0;
}

ラスト。main関数です。
とりあえずいつも通り画像を読み込んで二値化していきます。

その後処理の継続判定用の変数としてflagを定義します。
whileループを作りその中で処理αと処理βを行います。
その際に返ってくる値の和をflagに代入します。そのflagの値を評価して0以上ならば継続（どちらかの処理で変更があった）、0ならば停止します（どちらの処理でも変更がなかった）。

最後に細線化された画像を書き込んで終了です！長かった！

実際に使ってみた

コードを回してみました。
main関数にあるとおり、更新の都度の状態を表示するようにしたのでその経過を貼っていきますー

ここで処理が終了しました。
うーーん、、精度は悪いですね。細線化はできているものの細すぎてところどころ途切れていたりなくなっていたりしますね、、
まあそもそもこのZhang-Suenアルゴリズムも完璧ではないっぽくて割と癖があるそうなので仕方ないのかも、、

とりあえず今回はコレでおしまいです。
細線化は他にもアルゴリズムがあるのでまた記事にしようかと思います。

全体のコード

#include <iostream>
#include <opencv4/opencv2/opencv.hpp>
#define WIDTH 640
#define HEIGHT 480
using namespace std;


enum p_num{
    p2,p3,p4,p5,p6,p7,p8,p9
};

void readImage( const char *fname, unsigned char img[HEIGHT][WIDTH] )
{
    cv::Mat ReadImage( HEIGHT, WIDTH, CV_8U );

    ReadImage = cv::imread( fname,0);

    for( int w=0; w<WIDTH; w++ ){
        for( int h=0; h<HEIGHT; h++ )
        {
            img[h][w] = ReadImage.data[w+h*WIDTH];
        }
    }
}

void showImage(unsigned char img[HEIGHT][WIDTH]){
    cv::Mat saveImage( HEIGHT, WIDTH, CV_8U );

    for( int w=0; w<WIDTH; w++ ){
        for( int h=0; h<HEIGHT; h++ )
        {
            saveImage.data[w+h*WIDTH] = img[h][w];
        }
    }
    
    cv::imshow("aaa", saveImage);
    cv::waitKey(0);
}

void writeImage( const char *fname, unsigned char img[HEIGHT][WIDTH])
{
    cv::Mat saveImage( HEIGHT, WIDTH, CV_8U );

    for( int w=0; w<WIDTH; w++ ){
        for( int h=0; h<HEIGHT; h++ )
        {
            saveImage.data[w+h*WIDTH] = img[h][w];
        }
    }

    cv::imwrite( fname , saveImage );
}

void binary(unsigned char org[HEIGHT][WIDTH], unsigned char byn[HEIGHT][WIDTH]){
    for(int w=0; w<WIDTH; w++){
        for(int h=0; h<HEIGHT; h++){
            if((int)org[h][w] < 200)
                byn[h][w] = 0;
            else
                byn[h][w] = 255;
        }
    }
}

void f_cre(int y, int x, unsigned char img[HEIGHT][WIDTH], int f[9]){
    //この関数でchar型がint型にキャストされる
    f[0] = img[y - 1][x];
    f[1] = img[y - 1][x + 1];
    f[2] = img[y][x + 1];
    f[3] = img[y + 1][x + 1];
    f[4] = img[y + 1][x];
    f[5] = img[y + 1][x - 1];
    f[6] = img[y][x - 1];
    f[7] = img[y - 1][x - 1];
    f[8] = img[y - 1][x];
    
}

int A_P(int f[9]){
    int count = 0;
    for (int i = 0; i < 8; i++) {
        if ((f[i+1]-f[i]) == -255)
            count++;
    }
    //cout<<"A:"<<count<<endl;
    return count;
}

int B_P(int f[9]){
    int count = 0;
    for (int i = 0; i < 8; i++) {
        if (f[i] == 0)
            count++;
    }
    //cout<<"B:"<<count<<endl;
    return count;
}

bool judge_alpha(int y,int x,unsigned char img[HEIGHT][WIDTH]){
    int f[9]={0};
    f_cre(y, x, img, f);
    int A = A_P(f);
    int B = B_P(f);

    if((A==1)&&(2<=B)&&(B<=6)&&((f[p2]*f[p4]*f[p6])==0)&&((f[p4]*f[p6]*f[p8])==0))
        return true;
    
    return false;
}

bool judge_beta(int y,int x,unsigned char img[HEIGHT][WIDTH]){
    int f[9]={0};
    f_cre(y, x, img, f);
    int A = A_P(f);
    int B = B_P(f);

    if((A==1)&&(2<=B)&&(B<=6)&&((f[p2]*f[p4]*f[p8])==0)&&((f[p2]*f[p6]*f[p8])==0))
        return true;
    
    return false;
}

void change(unsigned char img[HEIGHT][WIDTH], bool reco[HEIGHT][WIDTH]) {
    for (int w = 1; w < (WIDTH - 1); w++){
        for (int h = 1; h < (HEIGHT - 1); h++){
            if (reco[h][w])
                img[h][w] = 255;
        }
    }
}


int alpha(unsigned char img[HEIGHT][WIDTH]) {
    bool reco[HEIGHT][WIDTH] = { false };
    int flag = 0;

    for (int h = 1; h < (HEIGHT - 1); h++){
        for (int w = 1; w < (WIDTH - 1); w++) {
            if((int)img[h][w]==0){
                if(judge_alpha(h, w, img)){
                    flag++;//画素の変更が起きた時にインクリメント
                    reco[h][w] = true;//画素の変更が起きる場所を記録
                }
            }
        }
    }
    change(img, reco);
    return flag;
}

int beta(unsigned char img[HEIGHT][WIDTH]) {
    bool reco[HEIGHT][WIDTH] = { false };
    int flag = 0;

    for (int h = 1; h < (HEIGHT - 1); h++){
        for (int w = 1; w < (WIDTH - 1); w++) {
            if((int)img[h][w]==0){
                if(judge_beta(h, w, img)){
                    flag++;//画素の変更が起きた時にインクリメント
                    reco[h][w] = true;//画素の変更が起きる場所を記録
                }
            }
        }
    }
    change(img, reco);
    return flag;
}
int main(int argc, const char * argv[]) {
    unsigned char img[HEIGHT][WIDTH]={0};
    cv::Mat saveImage( HEIGHT, WIDTH, CV_8U );
    readImage("thinning_5.png", img);
    
    binary(img, img);
    writeImage("gray_thin.png", img);

    
    int flag = 1;
    while (flag) {
        flag = 0;
        int a = alpha(img);//変更があったら1以上になる
        int b = beta(img);//変更があったら1以上になる
        flag = a + b;//どちらにも変更がない場合のみ０
        cout<<flag<<endl;
        showImage(img);
    }
    
    writeImage("thinning.png", img);
    return 0;
}