2D プロット・グラフのデータ・マークの数値化 (その4)

目盛の数値を読むことができたので、 データ・マークの位置を目盛に合わせて読み取れるようにします。

tociyuki/unplot2d-cxx11: from plot2d to numerical values (experimental works)

その 3 でオングストロームを認識できてなかったので、 横軸をナノメータで書き直したグラフから読み取っています。

読み取った実数値をカンマ区切り形式で出力しています。

"d/nm","Tc/K"
0.991561,2.49928
1.00211,2.56151
1.49789,2.46165
1.99367,2.12446
2.5,1.99855
3.99789,1.97395
5,2.01158
10,1.99855

その 1 から 3 までを一つにまとめます。 その 1 で試したデータ・マークの読み取りをメソッドに追加し、 さらにその 2 で試した X 軸の目盛位置の読み取りをメソッドにします。 X 軸と同じように Y 軸の目盛の読み取りもメソッドにします。

struct plot2d_type {
    // 略
    cv::Mat mark;
    cv::vector<cv::Point> data_mark;
    std::vector<int> xaxis_tick;
    std::vector<int> yaxis_tick;
    double xaxis_grad;
    double xaxis_offset;
    double yaxis_grad;
    double yaxis_offset;
    std::vector<cv::Point2d> data;

    // 略
    void choose_mark (cv::Mat& img);
    void match_data_mark (cv::Mat& img);

    bool detect_xaxis_ticks (cv::Mat& bin_img);
    bool detect_yaxis_ticks (cv::Mat& bin_img);
    void score_xaxis_ticks (cv::Mat& bin_img,
        std::vector<int>& peaks, std::vector<int>& peak_idx);
    void score_yaxis_ticks (cv::Mat& bin_img,
        std::vector<int>& peaks, std::vector<int>& peak_idx);

    double average_axis_interval (int frame_size,
        std::vector<int>& peaks, std::vector<int>& peak_idx);
    int min_axis_interval (int interval_min,
        std::vector<int>& peaks, std::vector<int>& peak_idx);
    void find_peaks (std::vector<int>& data, std::vector<int>& peaks);
    double correlation (cv::Mat& mat0, cv::Mat& mat1);

    void major_xaxis (void);
    void major_yaxis (void);
    void major_data_mark (void);
    bool at_peak (cv::Mat& img, int x, int y, float threshold);
};

int
main (int argc, char* argv[])
{
    plot2d_type plot2d;
    cv::Mat gray_img;
    cv::Mat bin_img;

    if (argc != 2) {
        std::cerr << "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;

    plot2d.detect_frame (bin_img);
    if (plot2d.frame.size () != 2) {
        std::cerr << "cannot detect plot frame" << std::endl;
        return EXIT_FAILURE;
    }

    plot2d.choose_mark (bin_img);
    plot2d.match_data_mark (bin_img);

    plot2d.detect_xaxis_ticks (bin_img);
    plot2d.detect_yaxis_ticks (bin_img);

    plot2d.detect_yaxis_label_rect (gray_img);
    plot2d.detect_xaxis_label_rect (gray_img);
    plot2d.detect_scale_label_rect (gray_img);

    tesseract::TessBaseAPI tess;
    plot2d.read_yaxis_label (tess, gray_img);
    plot2d.read_xaxis_label (tess, gray_img);
    plot2d.read_scale_label (tess, gray_img);

    plot2d.major_xaxis ();
    plot2d.major_yaxis ();
    plot2d.major_data_mark ();

    cv::Mat dst = img.clone ();
    cv::rectangle (dst, plot2d.frame[1], cv::Scalar (0, 0, 255), 1, 8);
    for (int x : plot2d.xaxis_tick) {
        cv::line (dst, cv::Point (x - 1, plot2d.frame[1].br ().y - 48),
                       cv::Point (x - 1, plot2d.frame[1].br ().y - 32),
                       cv::Scalar (0, 0, 255), 2, 8);
    }
    for (int y : plot2d.yaxis_tick) {
        cv::line (dst, cv::Point (plot2d.frame[1].x + 32, y - 1),
                       cv::Point (plot2d.frame[1].x + 48, y - 1),
                       cv::Scalar (0, 0, 255), 2, 8);
    }
    if (! plot2d.yaxis_label.empty ()) {
        cv::rectangle (dst, plot2d.yaxis_label_rect, cv::Scalar (0, 0, 255), 1, 8);
    }
    if (! plot2d.xaxis_label.empty ()) {
        cv::rectangle (dst, plot2d.xaxis_label_rect, cv::Scalar (0, 0, 255), 1, 8);
    }
    for (int i = 0; i < plot2d.scale_label.size (); ++i) {
        cv::rectangle (dst, plot2d.scale_label_rect[i], cv::Scalar (0, 0, 255), 1, 8);
    }
    for (int i = 0; i < plot2d.data_mark.size (); ++i) {
        cv::Rect r (plot2d.data_mark[i].x - plot2d.mark.cols / 2,
                    plot2d.data_mark[i].y - plot2d.mark.rows / 2,
                    plot2d.mark.cols, plot2d.mark.rows);
        cv::rectangle (dst, r, cv::Scalar (0, 255, 0), 1, 8);
    }

    std::cout << "\"" << plot2d.xaxis_label << "\","
              << "\"" << plot2d.yaxis_label << "\"" << std::endl;
    for (int i = 0; i < plot2d.data.size (); ++i) {
        std::cout << plot2d.data[i].x << "," << plot2d.data[i].y << std::endl;
    }
    //cv::imwrite ("result.png", dst);
    cv::namedWindow ("PLOT", cv::WINDOW_NORMAL);
    cv::imshow ("PLOT", dst);
    cv::waitKey ();

    return EXIT_SUCCESS;
}

データ・マークの読み取りは、 choose_markmatch_data_mark の 2 つのメソッドに分けます。 前者は枠内の輪郭からマークを一つ選びます。 後者は選んだマークでテンプレート・マッチングをして data_mark メンバに画像の座標系で位置を記録します。 その 1 からの改良点として、 data_mark にはマークの中心位置を格納するように変更しました。 at_peak は関数をメソッドにしただけで中身は前と同じです。

void
plot2d_type::choose_mark (cv::Mat& img)
{
    for (int i = 0; i < bound.size (); ++i) {
        if (frame[0] == bound[i]
                || (frame[0] & bound[i]) != bound[i]
                || abs (bound[i].width - bound[i].height) > 8
                || std::max (bound[i].width, bound[i].height) > 36)
            continue;
        mark = img (bound[i]).clone ();
        break;
    }
}

void
plot2d_type::match_data_mark (cv::Mat& img)
{
    int const xoffset = mark.cols / 2;
    int const yoffset = mark.rows / 2;
    cv::Mat result;
    cv::matchTemplate (img, mark, result, cv::TM_CCORR_NORMED);
    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)) {
                data_mark.push_back (cv::Point (x + xoffset, y + yoffset));
            }
        }
    }
}

X 軸の目盛位置の読み取りは、 その 2 と同じで引数のいくつかをメンバにしているだけなので、 省略します。 Y 軸の目盛位置の読み取りは、 X 軸同様に、 ヒストグラムを作り、 それのピーク位置を find_peaks メソッドでリストアップします。 その中から目盛らしきものを抽出して、 均等に目盛を軸に割り当てます。 find_peaks メソッド、 average_axis_interval メソッド、 correlation メソッドの 3 つは、 単に関数をメソッドにしただけなので省略します。

bool
plot2d_type::detect_yaxis_ticks (cv::Mat& bin_img)
{
    cv::Rect axis_rect (frame[0].x, frame[0].y,
        frame[1].x + 24 - frame[0].x, frame[0].height);
    cv::Mat axis = bin_img (axis_rect);
    std::vector<int> histogram (axis.rows, 0);
    for (int i = 0; i < axis.rows; ++i) {
        histogram[i] = cv::countNonZero (axis (cv::Rect (0, i, axis.cols, 1)));
    }
    std::vector<int> peaks;
    find_peaks (histogram, peaks);
    std::vector<int> peak_idx;
    score_yaxis_ticks (bin_img, peaks, peak_idx);
    double const interval = average_axis_interval (frame[1].height, peaks, peak_idx);
    double const y0 = peaks[peak_idx[0]] + frame[0].y;
    int const i0 = (int)(((double)frame[1].y - y0) / interval);
    for (int i = i0; y0 + interval * i <= (double)frame[1].br ().y + 8; ++i) {
        int pos = cvRound (y0 + interval * i);
        if (pos >= frame[1].y) {
            yaxis_tick.push_back (pos);
        }
    }
    return true;
}

void
plot2d_type::score_yaxis_ticks (cv::Mat& bin_img,
                    std::vector<int>& peaks,
                    std::vector<int>& peak_idx)
{
    for (int i = 0; i < peaks.size (); ++i) {
        cv::Rect tr (frame[0].x, frame[0].y + peaks[i] - 18,
                     frame[1].x + 24 - frame[0].x, 36);                     
        cv::Mat tpl = bin_img (tr);
        int score = 0;
        for (int j = 0; j < peaks.size (); ++j) {
            if (i == j) continue;
            cv::Rect ir (frame[0].x, frame[0].y + peaks[j] - 18, tr.width, tr.height);
            cv::Mat roi = bin_img (ir);
            double corr = correlation (tpl, roi);
            if (corr > 0.92) {
                ++score;
            }
        }
        if (score * 3 > peaks.size ()) {
            peak_idx.push_back (i);
        }
    }
}

画像上の目盛位置を目盛の数字に対応付けて、 画像上の位置から軸で指定した値へ変換するパラメータを求めます。 今はリニア・スケールの軸だけを扱っているので、 X 軸なら a * x + b の式の ab を線形最小二乗法で求めます。

struct linear_least_square_method_type {
    int n;
    double sum_x;
    double sum_xx;
    double sum_y;
    double sum_yx;
    double grad;
    double offset;
    void clear (void);
    void put (double x, double y);
    void finalize (void);
};

X 軸の目盛に対して変換式のパラメータを計算します。 X 軸の目盛の画像位置と、 目盛の数字列の横方向の中心が近いときに、 目盛と数値が対応しているとみなし、 積算を進めます。 積算が終わったら、 傾きとオフセットを求めてメンバに記録します。

void
plot2d_type::major_xaxis (void)
{
    linear_least_square_method_type lsq;
    lsq.clear ();
    int const tick_w = (xaxis_tick[1] - xaxis_tick[0]) / 2;
    for (int i = 0; i < xaxis_tick.size (); ++i) {
        int const x = xaxis_tick[i];
        for (int j = 0; j < scale_label_rect.size (); ++j) {
            int const scale_w = scale_label_rect[j].width / 2;
            int const scale_x = scale_label_rect[j].x + scale_w;
            int const scale_y = scale_label_rect[j].y;
            if (abs (scale_x - x) < tick_w && frame[0].br ().y < scale_y) {
                lsq.put (x, std::stod (scale_label[j]));
            }
        }
    }
    lsq.finalize ();
    xaxis_grad = lsq.grad;
    xaxis_offset = lsq.offset;
}

Y 軸も同様に目盛と数字列の対応を付けて、 傾きとオフセットを求めてメンバに記録します。 やはり、 Y 軸の目盛の画像位置と、 目盛の数字列の縦方向の中心が近いときに両方が近いと判断しています。

void
plot2d_type::major_yaxis (void)
{
    linear_least_square_method_type lsq;
    lsq.clear ();
    int const tick_h = (yaxis_tick[1] - yaxis_tick[0]) / 2;
    for (int i = 0; i < yaxis_tick.size (); ++i) {
        int const y = yaxis_tick[i];
        for (int j = 0; j < scale_label_rect.size (); ++j) {
            int const scale_h = scale_label_rect[j].height / 2;
            int const scale_x = scale_label_rect[j].br ().x;
            int const scale_y = scale_label_rect[j].y + scale_h;
            if (abs (scale_y - y) < tick_h && scale_x < frame[0].x) {
                lsq.put (y, std::stod (scale_label[j]));
            }
        }
    }
    lsq.finalize ();
    yaxis_grad = lsq.grad;
    yaxis_offset = lsq.offset;
}

両軸の変換式の係数が求まったので、 データ・マークの画像座標位置から実数値へ変換して data メンバへ記録します。

void
plot2d_type::major_data_mark (void)
{
    for (int i = 0; i < data_mark.size (); ++i) {
        double const x = xaxis_grad * data_mark[i].x + xaxis_offset;
        double const y = yaxis_grad * data_mark[i].y + yaxis_offset;
        data.push_back (cv::Point2d (x, y));
    }
    std::sort (data.begin (), data.end (), [](cv::Point2d& a, cv::Point2d& b) { return a.x < b.x; });
}

線形最小二乗法の定義通りに 2 つの係数の推定値を求めます。

void
linear_least_square_method_type::clear (void)
{
    n = 0;
    sum_x  = 0.0;
    sum_xx = 0.0;
    sum_y  = 0.0;
    sum_yx = 0.0;
}

void
linear_least_square_method_type::put (double x, double y)
{
    ++n;
    sum_x  += x;
    sum_xx += x * x;
    sum_y  += y;
    sum_yx += y * x;
}

void
linear_least_square_method_type::finalize (void)
{
    double const d = sum_xx * n - sum_x * sum_x;
    grad   = (sum_yx * n - sum_x * sum_y) / d;
    offset = (sum_xx * sum_y - sum_yx * sum_x) / d;
}