2D プロット・グラフの枠線の認識 (その1)

これは、 主に科学技術論文に掲載されているグラフから数値を復元する試みの初めの一歩です。 初めの一歩だけに、 問題をシンプルにして、 PNG に描いたモノクロ 2 値ロスレス圧縮・歪みなし・モアレなし・ノイズなしという最も扱いやすいグラフの画像からの数値読み取りを試みます。

また、 画像の解像度が 300 dpi であるとして、 ピクセル数をハードコードしています。

今回は、 グラフの枠線の認識と、 その中に点在するマークを画像処理のテンプレート・マッチングで認識し、 枠線の大きさでゼロから 1 の間に正規化した数値を求めるところまで試します。 しかも、 マークは 1 種類だけとして、 マーク同士の重なりもなし、 といった具合に認識に苦労することがない画像を扱います。 正規化した読み出しだけなので、 枠線とマークを OpenCV 2.4 で認識させます。 そのうち、 文字認識が必要になるので、 近いうちに OpenCV 3 へ移行する予定ですが、 まずは使い慣れている 2.4 で認識をおこなっていきます。

都合の良いことに、 多くの科学技術論文に掲載されているグラフのほとんどは、 長方形の枠線でプロット・エリアを囲んであります。 従って、 まずやるべきことは枠線を認識することです。 ただし、 枠線は単なる長方形ではなく、 目盛が重なっていたり、 データのマークが重なっていることもあります。 それを考慮して、 枠線を認識するため、 OpenCVcv::findContour で輪郭線抽出をおこなっておき、 その中から、 画像の 4 分の 1 以上をしめる大きさで、 なおかつ、 cv::approxPolyDP で多角形近似をして四角形になる輪郭を探しだす方針でいきます。

作るのはコマンドライン・アプリケーションとし、 コマンドライン引数に画像ファイル名を与えることにします。 それを読み込んで、 グレー・スケールにし、 2 値化画像に変換していく手順は、 定番通りです。 コンパイルして枠とマークを認識させてみます。

$ clang++ -std=c++11 `pkg-config --cflags opencv` frame.cpp -o frame `pkg-config --libs opencv` -lm
$ ./frame plot.png
0.100105 0.550289
0.100105 0.500578
0.149631 0.47052
0.200211 0.201156
0.499473 0.110983
0.249737 0.100578
0.998946 0.100578
0.399368 0.0809249
$

認識した枠は赤線で、 マークは緑の長方形で囲んで白黒反転したグラフ画像に重ね書きしたものを下に示します。 読み取りに使ったグラフは、 私の修士論文からで、 金属バナジウムと非晶質シリコンの多層膜の超伝導転移温度のグラフです。 バナジウム層の厚さは 30 Å、 グラフの横軸は非晶質シリコン層の厚さで、 縦軸は超伝導転移温度です。 成膜方法は超高真空蒸着法。 4 端子法によるシート電気伝導率測定。 温度制御はキニーポンプによるヘリウム減圧。

今後、 目盛の抽出、 目盛の数字列読み取りができるようになると、 グラフ上の数値へ変換することができるようになるのでしょう。 さらにラベルを読めるようになれば、 単位を認識できることになります。

#include "opencv/cv.h"
#include "opencv/highgui.h"
#include <cstdlib>
#include <iostream>
#include <algorithm>

void detect_frame (cv::Mat& img,
              cv::vector<cv::vector<cv::Point> >& contours,
              int i,
              cv::vector<cv::Vec4i>& hierarchy,
              cv::vector<cv::Rect>& frame);

bool rect_contains (cv::Rect const& a, cv::Rect const& b);

bool at_peak (cv::Mat& result, int x, int y, float threshold);


int
main (int argc, char* argv[])
{
    cv::Mat gray_img;
    cv::Mat bin_img;
    cv::Mat dst;
    cv::vector<cv::vector<cv::Point> > contours;
    cv::vector<cv::Vec4i> hierarchy;
    cv::vector<cv::Rect> bound;
    cv::vector<cv::Rect> frame;

    if (argc != 2) {
        std::cout << "usage: " << argv[0] << " image-file-name" << std::endl;
        return EXIT_FAILURE;
    }
    cv::Mat img = cv::imread (argv[1], cv::IMREAD_UNCHANGED);
    if (! img.data)
        return EXIT_FAILURE;
    cv::cvtColor (img, gray_img, CV_BGR2GRAY);
    cv::threshold (gray_img, bin_img, 0, 255, cv::THRESH_BINARY|cv::THRESH_OTSU);
    bin_img = ~bin_img;
    cv::cvtColor (bin_img, dst, CV_GRAY2BGR);

//@<輪郭線を抽出する@>
//@<輪郭線から枠を探し出す@>
//@<枠の中からマークを一個無作為に選んでテンプレート画像にする@>
//@<テンプレート・マッチングをおこなう@>
//@<枠の中でテンプレート・マッチング結果のピークを探す@>
//@<ピーク位置を枠中の正規化座標に変換して出力する@>

    cv::namedWindow ("IMAGE", cv::WINDOW_NORMAL);
    cv::imshow ("IMAGE", dst);
    cv::waitKey ();
    return EXIT_SUCCESS;
}

//@<枠の判定をおこなう detect_frame を定義する@>
//@<四角形の包含関係をチェックする rect_contains を定義する@>
//@<閾値以上のピーク位置かどうかチェックする at_peak を定義する@>

輪郭線抽出は cv::findContour でおこないます。 この手続きは画像を書き直してしまうので、 マークの認識に使うための 2 値化画像を複製しておきます。 輪郭線の階層構造を 2 段階で求める CV_RETR_CCOMP を指定しています。 これにより枠とマークそれぞれを親階層の輪郭線で求めることができ、 さらに枠の中抜きの内側の輪郭線を子階層で得ることができるようになります。

//@<輪郭線を抽出する@>=
    cv::Mat ctr_img = bin_img.clone ();
    cv::findContours (ctr_img, contours, hierarchy, CV_RETR_CCOMP, CV_CHAIN_APPROX_SIMPLE);

マーク、 枠、 目盛の文字、 ラベルの文字のそれぞれの輪郭線が抽出されるので、 それぞれの外接長方形を記録しておきます。 さらに、 抽出された中から枠の輪郭線を探します。

//@<輪郭線から枠を探し出す@>=
    for (int i = 0; i >= 0; i = hierarchy[i][0]) {
        bound.push_back (cv::boundingRect (contours[i]));
        detect_frame (img, contours, i, hierarchy, frame);
    }
    if (frame.size () != 2) {
        std::cerr << "cannot detect plot frame" << std::endl;
        return EXIT_FAILURE;
    }

輪郭線ベクタ contoursi 要素が、 枠かどうかを判定します。 グラフの枠は画像の 4 分の 1 より大きいはずなので、 面積を cv::contourArea で求めます。 続いて cv::approxPolyDP で多角形近似をおこない、 目盛や重ね書くされたマークといった小さな付属図形を取り除きます。 その結果、 頂点数が 4 つになると枠線の可能性が高くなります。 さらに、 枠は線分の長方形であるため、 輪郭線の階層構造で子要素をもっているはずなので、 それがあるかどうかを調べます。 子階層の中から、 やはり大きな面積の輪郭線を探し出します。 内側の輪郭線にも目盛とマークが重ね書きされていることがあるので、 これも多角形近似して四角形になることを確かめます。 なお、 この後に、 枠の線幅をチェックした方が良いのでしょうが、 ここでは省いています。

frame ベクタには、 先頭要素に枠線とそれに重ね書きされたマークや目盛を含んだ外接長方形が入ります。 この外接長方形の中にマークが存在しているはずなので、 テンプレート・マチング結果のピークを、 この中で探すようにします。 もう一つ、 frame ベクタに、 四角形近似した枠の内側と外側の中心を通る四角形を格納しておきます。 こちらはマークの座標の正規化に使います。

なお、 この枠の判定は手抜きで、 四角形が長方形かどうか調べていません。 今後の課題で、 頂点の角度が π/2 かどうかの判定をおこなうように改善するべきでしょう。 (グラフに理論曲線が含まれているとき、 下のやりかたでは枠の内線が正しく求まらない場合があることに気がついたので、 その3 で計算方法を変更しています)

//@<枠の判定をおこなう detect_frame を定義する@>=
void
detect_frame (cv::Mat& img,
              cv::vector<cv::vector<cv::Point> >& contours,
              int const i,
              cv::vector<cv::Vec4i>& hierarchy,
              cv::vector<cv::Rect>& frame)
{
    double area = cv::contourArea (contours[i], false);
    if (area <= img.cols * img.rows / 4)
        return;
    cv::vector<cv::Point> approx;
    double epsilon = 0.01 * cv::arcLength (contours[i], true);
    cv::approxPolyDP (cv::Mat (contours[i]), approx, epsilon, true);
    if (approx.size () != 4)
        return;
    if (hierarchy[i][2] < 0)
        return;
    for (int i1 = hierarchy[i][2]; i1 >= 0; i1 = hierarchy[i1][0]) {
        double area1 = cv::contourArea (contours[i1], false);
        if (area1 <= img.cols * img.rows / 4)
            continue;
        cv::vector<cv::Point> approx1;
        double epsilon1 = 0.01 * cv::arcLength (contours[i1], true);
        cv::approxPolyDP (cv::Mat (contours[i1]), approx1, epsilon1, true);
        if (approx1.size () == 4) {
            frame.push_back (cv::boundingRect (contours[i]));
            cv::Rect r0 = cv::boundingRect (approx);
            cv::Rect r1 = cv::boundingRect (approx1);
            cv::Rect r ((r0.x + r1.x) / 2, (r0.y + r1.y) / 2,
                        (r0.width + r1.width) / 2,
                        (r0.height + r1.height) / 2);
            frame.push_back (r);
            return;
        }
    }
    return;
}

枠が求まり、 グラフのプロット領域がわかったので、 その中のマークを抽出します。 マークは共通パターンで繰り返し描画されているはずであり、 回転も変形もないはずなので、 テンプレート・マッチングで特徴抽出してみることにします。 そのためには、 マークのテンプレートが必要なので、 枠の中にある小さな輪郭線を一つ選びます。 今回は、 マークは 1 種類しかないと仮定しているので、 無作為に一つ選んでいます。 その際、 これまで読んだ論文で 3 mm を越える大きさのマークを使っているグラフはなかった覚えがあるため、 300dpi での 3 mm、 すなわち 36 ピクセル以下で幅と高さが同じ程度の輪郭線を選びます。

//@<枠の中からマークを一個無作為に選んでテンプレート画像にする@>=
    cv::Mat mark;
    for (int i = 0; i < bound.size (); ++i) {
        // 36 px == 3mm * 300 px / 25.4 mm
        if (! rect_contains (frame[0], bound[i])
                || abs (bound[i].width - bound[i].height) > 8
                || std::max (bound[i].width, bound[i].height) > 36)
            continue;
        mark = bin_img (bound[i]).clone ();
    }

上で使う、 四角形 a の中に四角形 b が内包されているときに真を返す関数を作っておきます。

//@<四角形の包含関係をチェックする rect_contains を定義する@>=
bool
rect_contains (cv::Rect const& a, cv::Rect const& b)
{
    return a.x < b.x && a.y < b.y
        && b.br ().x < a.br ().x && b.br ().y < a.br ().y;
}

選んだマークを使ってテンプレート・マッチングをおこないます。 マッチングの結果は、 ピクセルごとに、 そのピクセルがテンプレートに右上になる尤もらしさを正規化した 2 次元マップになります。 マッチング方式に cv::TM_CCORR_NORMED を選んだときは、 ゼロから 1.0 の間で 1.0 に近いほど尤もらしいことになります。 枠の中で、 しきい値以上の値で 2 次元のピークになっているピクセル座標を選び出します。

//@<枠の中でテンプレート・マッチング結果のピークを探す@>=
    cv::Mat result;
    cv::matchTemplate (bin_img, mark, result, cv::TM_CCORR_NORMED);
    cv::vector<cv::Point> plots;
    float threshold = 0.7f;
    for (int y = 0; y < result.rows; ++y) {
        for (int x = 0; x < result.cols; ++x) {
            if (frame[0].contains (cv::Point (x, y))
                    && at_peak (result, x, y, threshold)) {
                plots.push_back (cv::Point (x, y));
            }
        }
    }

あるピクセル座標が閾値以上のピークの座標に一致しているかどうかを判定します。 今回は、 ピークが複数のピクセルにまたがっているときは、 左上の座標を選ぶようにしています。 それよりも、 意味合いからして重心位置にあるかどうかを判定する方が良いかもしれません。

//@<閾値以上のピーク位置かどうかチェックする at_peak を定義する@>=
bool
at_peak (cv::Mat& result, int x, int y, float threshold)
{
    if (result.at<float> (y, x) < threshold)
        return false;
    float m = 0.0f;
    int xp = -1, yp = -1;
    for (int dy = -3; dy < 4; ++dy) {
        for (int dx = -3; dx < 4; ++dx) {
            if (result.at<float> (y + dy, x + dx) > m) {
                m = result.at<float> (y + dy, x + dx);
                yp = y + dy;
                xp = x + dx;
            }
        }
    }
    return xp == x && yp == y;
}

これで、 枠とマークのピクセル座標が求まったので、 マークの位置を正規化してリストアップします。

//@<ピーク位置を枠中の正規化座標に変換して出力する@>=
    cv::rectangle (dst, frame[1].tl (), frame[1].br (), cv::Scalar (0, 0, 255), 2, 8);
    float mark_mag = mark.cols / 2;
    for (int i = 0; i < plots.size (); ++i) {
        cv::Point p1 = plots[i];
        cv::Point p2 (p1.x + mark.cols, p1.y + mark.rows); 
        cv::rectangle (dst, p1, p2, cv::Scalar (0, 255, 0), 2, 8);

        float x = ((float)p1.x + mark_mag - (float)frame[1].x) / (float)frame[1].width;
        float y = ((float)frame[1].br ().y - (float)p1.y - mark_mag) / (float)frame[1].height;
        std::cout << x << " " << y << std::endl;
    }

    cv::namedWindow ("IMAGE", cv::WINDOW_NORMAL);
    cv::imshow ("IMAGE", dst);
    cv::waitKey ();
    return EXIT_SUCCESS;
}

これは、あくまでも理想的な場合の第一歩にすぎず、 考慮しなければならないことが山ほどあります。 今回はマークが 1 種類しかないと仮定していましたが、 多くのグラフでは数種類のマークが混在するので、 判別できるようにしなければなりません。 さらに、 マーク同士が重なることがあるので、 それも分離して抽出できなければなりません。 マークの間に理論式を当てはめた曲線が引いてあることも良くあることで、 曲線とマークの分離もできるようにならないといけません。 マークにエラーバーが付いているときもあるので、 それも認識できないといけません。 グラフ中に、 矢印やら注釈文字列が混じっていることもあるので、 それらを分離する方法が必要です。 中にはグラフ中に説明のための図も入っていることがあるので、 それらをマークとして誤認識しないようにしないといけません。

そこまでやってもまだ途上です。 理想的な画像ではなく、 PDF 埋め込みの JPEG 画像からグラフであることを認識して、 数値化処理をするかどうか判定できるようにしないといけません。 それだけで済めばまだ良いわけで、 印刷雑誌からスキャンした、 ノイズ混じり歪みや回転された画像からグラフを読むところまで辿り着けて、 ようやく実用の手前に立つことができます。