2D プロット・グラフの横軸目盛の認識 (その2)

グラフの枠の付属物である目盛の位置を検出します。 枠の下辺である横軸がリニア・スケールであるとして、 目盛を読んでみます。

画像の特徴としては、 等間隔で同じ形状の目盛が並ぶこと、 目盛の個数が少ないこと、 目盛が短いこと、 データ・マークが大きなノイズになることがあげることができます。 等間隔で並ぶことを利用しようと、 当初、考えた方法は、 枠の下枠近傍の狭い領域でヒストグラムを作り、 それの自己相関関数を求めるやりかただったのですが、 枠に重なっているデータのマークは系統誤差であるため、 自己相関関数を求めても目盛位置のピークだけが残るようなことはなく、 この案は没になりました。

今回の手順です。

  1. ヒストグラムを作ります。
  2. ヒストグラムのピーク位置を検出します。
  3. ヒストグラムのピーク近傍の小さな画像をテンプレートにして、 他のピーク位置との相関関数を計算し、 スコアをつけます。
  4. スコアから、 目盛の可能性が高いピーク位置を抜き出します。
  5. ピーク位置から、 最小の等間隔を求めます。
  6. 最小の等間隔で、 目盛位置を予測します。

目盛の形状の特徴を利用せずに位置を検出しており、 とりあえず読めたというだけのものです。 これも近いうちに没にして、 さらに良い方法を見つけるつもりです。

前回のやりかたで求めた枠の四角形ベクタを使って、 X 軸の目盛の位置を検出します。 main 関数の枠の検出直後に次のコード片を挿入します。

    std::vector<int> xticks_pos;
    std::vector<int> xticks_pos;
    detect_xaxis_ticks (bin_img, frame, xticks_pos);

    for (int x : xticks_pos) {
        cv::line (dst, cv::Point (x - 1, frame[1].br ().y - 48),
                       cv::Point (x - 1, frame[1].br ().y - 32),
                       cv::Scalar (0, 0, 255), 2, 8);
    } 

枠の下辺近傍の細長い長方形領域の中で、 縦方向にゼロでないピクセルの個数を数えて、 ヒストグラムを作ります。 目盛の長さは数 mm が大半なので、 長方形の高さは内側に 2mm にしています。 外側は枠の輪郭線に外接する長方形の下辺にします。

bool
detect_xaxis_ticks (cv::Mat& bin_img,
                    cv::vector<cv::Rect>& frame,
                    std::vector<int>& ticks_pos)
{
    // x axis and ticks may be drawn near frame line
    // 24px == 2mm * 300px / 25.4mm
    cv::Rect axis_rect (frame[0].x, frame[1].br ().y - 24.0,
        frame[0].width, frame[0].br ().y - frame[1].br ().y + 24.0);
    cv::Mat axis = bin_img (axis_rect);
    std::vector<int> histogram (axis.cols, 0);
    for (int i = 0; i < axis.cols; ++i) {
        histogram[i] = cv::countNonZero (axis (cv::Rect (i, 0, 1, axis.rows)));
    }
    std::vector<int> peaks;
    find_peaks (histogram, peaks);
    std::vector<int> tick_peak_idx;
    score_xaxis_ticks (bin_img, frame, peaks, tick_peak_idx);
    double const interval = average_axis_interval (frame[1].width, peaks, tick_peak_idx);
    double const x0 = peaks[tick_peak_idx[0]] + frame[0].x;  // 2017-06-08 修正
    //double const x0 = peaks[tick_peak_idx[0]];
    int const i0 = (int)(((double)frame[1].x - x0) / interval);
    for (int i = i0; x0 + interval * i <= (double)frame[1].br ().x + 8; ++i) {
        int pos = cvRound (x0 + interval * i);
        if (pos >= frame[1].x) {
            ticks_pos.push_back (pos);
        }
    }
    return true;
}

ヒストグラムのピーク位置を求めます。 オーソドックスに Savitzky と Golay の平滑化微分法を使います。 左から右へと順に微分値を追っていき、 微分値がプラスからマイナスへ変化する場所がピーク位置の候補です。 途中にゼロが続くときがあるときは、 ゼロの中央をピーク位置の候補に選びます。 用心して、 候補位置の周辺で極大値の位置をチェックして、 そこをピーク位置として返します。 なお、 これの元は FORTRAN のコードだった名残りで、 データ配列の添字が 1 から始まる流儀になっています。 ベクタ要素を読む段階で 1 を引いて、 添字がゼロから始まる C 言語の流儀に合わせてあります。

void
find_peaks (std::vector<int>& data, std::vector<int>& peaks)
{
    int const q = 4;
    int sign = -1;
    int span = 0;
    for (int i = q + 1; i <= data.size () - q; ++i) {
        double p = 0.0;
        for (int j = 1; j <= q; ++j) {
            p += (double)(data[i + j - 1] - data[i - j - 1]) * j;
        }
        int pos = 0;
        int left_sign = sign;
        sign = p < 0 ? -1 : p > 0 ? +1 : 0;
        if (sign > 0) {
            span = 0;
        }
        else if (sign == 0) {
            ++span;
        }
        else if (left_sign == 0) {
            pos = i - 1 - span / 2;
        }
        else if (left_sign > 0) {
            pos = i;
        }
        else {
            span = 0;
        }
        if (left_sign >= 0 && sign < 0) {
            int peak_pos = pos;
            for (int k = pos - q; k <= pos + q; ++k) {
                if (data[k - 1] > data[peak_pos - 1]) {
                    peak_pos = k;
                }
            }
            if (data[peak_pos - 1] > 3) {
                peaks.push_back (peak_pos - 1);
            }
            i += q;
        }
    }
}

続いて、 枠や目盛にデータ・マークが重なっている場合をはじくために、 求めたピーク位置の周辺の小領域が目盛であるかどうかスコアで判断します。 目盛は同じパターンで繰り返し出現するはずなので、 他のピーク位置の同じ大きさの小領域との相関があるときにスコアをアップします。 そうして、 スコアがピークの位置の個数の 3 分の 1 を越えるときに、 目盛位置の候補とします。

void
score_xaxis_ticks (cv::Mat& bin_img,
                    cv::vector<cv::Rect>& frame,
                    std::vector<int>& peaks,
                    std::vector<int>& tick_peak_idx)
{
    for (int i = 0; i < peaks.size (); ++i) {
        cv::Rect tr (frame[0].x + peaks[i] - 18, frame[1].br ().y - 24.0,
                     36, frame[0].br ().y - frame[1].br ().y + 24.0);
        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 + peaks[j] - 18, frame[1].br ().y - 24.0,
                         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 ()) {
            tick_peak_idx.push_back (i);
        }
    }
}

小領域同士の相関は素朴なやりかたで書いています。

double
correlation (cv::Mat& mat0, cv::Mat& mat1)
{
    double den = 0.0, num0 = 0.0, num1 = 0.0;
    for (int dy = 0; dy < mat0.rows; ++dy)
        for (int dx = 0; dx < mat0.cols; ++dx) {
            double const z0 = mat0.at<unsigned char> (dy, dx);
            double const z1 = mat1.at<unsigned char> (dy, dx);
            den  += z0 * z1;
            num0 += z0 * z0;
            num1 += z1 * z1;
        }
    return den / std::sqrt (num0 * num1);
}

求めた目盛の候補位置から、 目盛の間隔を推定します。 最小間隔を求めておいて、 それの倍数の候補位置間隔を足していきます。 最小間隔換算の間隔の個数を同時に計数しておいて、 最小間隔の平均値を計算します。

double
average_axis_interval (int frame_size,
                        std::vector<int>& peaks,
                        std::vector<int>& tick_peak_idx)
{
    int interval_min = min_axis_interval (frame_size, peaks, tick_peak_idx);
    int interval_sum = 0;
    int interval_count = 0;
    for (int j = 1; j < tick_peak_idx.size (); ++j) {
        int const i0 = tick_peak_idx[j - 1];
        int const i1 = tick_peak_idx[j];
        int q = peaks[i1] - peaks[i0];
        if (q - interval_min < 2) {
            interval_sum += q;
            ++interval_count;
        }
        else {
            for (int i = 2; interval_min * i < frame_size; ++i) {
                if (abs (interval_min * i - q) < 2 * i) {
                    interval_sum += q;
                    interval_count += i;
                    break;
                }
            }
        }
    }
    return (double)interval_sum / (double)interval_count;
}

目盛の候補位置リストから、 最小間隔を求めます。

int
min_axis_interval (int interval_min,
                    std::vector<int>& peaks,
                    std::vector<int>& tick_peak_idx)
{
    for (int j = 1; j < tick_peak_idx.size (); ++j) {
        int const i0 = tick_peak_idx[j - 1];
        int const i1 = tick_peak_idx[j];
        int q = peaks[i1] - peaks[i0];
        if (q < interval_min) {
            interval_min = q;
        }
    }
    return interval_min;
}