読者です 読者をやめる 読者になる 読者になる

kivantium活動日記

プログラムを使っていろいろやります

画像処理で大腸菌のコロニー数を数えたかった話 その2

kivantium.hateblo.jp
の続編です。

前回は先行研究を一切調べずに適当にやったので今回は少し先行研究調査をして前回よりはまともな方法を実装したいと思います。

先行研究

Google検索したところ

ということが分かりました。

見た論文の中で一番簡単そうだったCell counting based on local intensity maxima grouping for in-situ microscopy - IEEE Xplore Documentを実装することにしました。

論文の内容と妥協ポイント

この論文では

  • Cell cluster Region Segmentation
  • Local Intensity Maxima Detection
  • Local Intensity Maxima Grouping
  • Cell Counting

の4段階で細胞数をカウントしていました。

Cell cluster Region SegmentationではMartinezの方法で細胞領域のセグメンテーションを行うと書いてあったのですが、Martinezの論文が入手できませんでした。そこでそれっぽいセグメンテーションができる閾値を手動で設定することにしました(妥協ポイント1)

Local Intensity Maxima Detectionではbilateral filterを掛けた後、各ピクセルについてそのピクセルを中心とする3x3の領域で最も明るい場合にlocal intensity maximumと判定するという処理をしていました。実験の結果3x3では小さすぎるようだったのでwindow sizeも手動で設定することにしました(妥協ポイント2)

Local Intensity Maxima Groupingではlocal intensity maximumを5つの基準に従ってグルーピングしていました。しかし、論文の基準は「2つのmaximumを結ぶ直線が細胞の境界と交差しない」のような実装が面倒なものが多かったので、2つのmaximumが同じwindowに入るならグループにするという単純な基準を採用しました(妥協ポイント3)

Cell Countingはグループの数を数える処理です。

実装

以上の妥協を行った実装が以下のようになります。(何も見直ししてないからコードが汚い)

#include <opencv2/core.hpp>
#include <opencv2/imgproc.hpp>
#include <opencv2/imgcodecs.hpp>
#include <opencv2/highgui.hpp>

#include <iostream>
#include <vector>
#include <map>
#include <cmath>
#include <algorithm>

cv::Mat org, src, bin, blured, dst;
int wsize = 5;
int threshold = 200;

using P = std::pair<int, int>;

struct UnionFind {
    std::vector<int> par;
    UnionFind(int n) : par(n, -1) {}
    int find(int x) { return par[x] < 0 ? x : par[x] = find(par[x]); }
    bool same(int x, int y) { return find(x) == find(y); }
    void unite(int x, int y) {
        x = find(x); y = find(y);
        if (x == y) return;
        if (par[y] < par[x]) std::swap(x, y);
        if (par[x] == par[y]) --par[x];
        par[y] = x;
    }
};

void drawWindow(void) {
    cv::Mat img = org.clone();
    cv::bilateralFilter(src, blured, 7, 5, 35);
    cv::threshold(blured, bin, threshold, 255, cv::THRESH_BINARY);
    std::vector<P> cell;
    for(int y=wsize; y<blured.cols-wsize; ++y) {
        for(int x=wsize; x<blured.rows-wsize; ++x) {
            int c = blured.data[y*blured.step+x]; 
            bool is_max = true;
            for(int dy=-wsize; dy<=wsize; ++dy) {
                for(int dx=-wsize; dx<=wsize; ++dx) {
                   if(blured.data[(y+dy)*blured.step+(x+dx)] > c) is_max = false;
                }
            }
            if(is_max && bin.data[y*bin.step+x] == 255) {
                int p = y*img.step+x*img.elemSize();
                img.data[p+0] = 0;
                img.data[p+1] = 0;
                img.data[p+2] = 255;
                cell.push_back(P(y, x));
            }
        }
    }

    UnionFind uf(cell.size());
    for(size_t i=0; i<cell.size(); ++i) {
        for(size_t j=0; j<cell.size(); ++j) {
            if(abs(cell[i].first - cell[j].first) < wsize && abs(cell[i].second - cell[j].second) < wsize) {
                uf.unite(i, j);
            }
        }
    }

    int pop = 0;
    for(size_t i=0; i<cell.size(); ++i) {
        if(uf.find(i) == i) pop++;
    }
    std::cout << "population is: " << pop << std::endl;

    cv::imshow("Source", img);
    cv::imshow("Binary", bin);

}

// トラックバーに変化があった時に呼ばれる関数
void on_trackbar( int, void* ) {
    drawWindow();
}

int main(int argc, const char** argv) {
    org = cv::imread(argv[1]);
    src = cv::imread(argv[1], cv::IMREAD_GRAYSCALE);
    dst = cv::Mat(src.size(), CV_8UC3);

    // 画像の読み込みに失敗したらエラー終了する
    if(src.empty()) {
        std::cerr << "Failed to open image file." << std::endl;
        return -1; 
    }

    // トラックバー付きウィンドウの作成
    cv::namedWindow("Binary");
    cv::createTrackbar("threshold", "Binary", &threshold, 255, on_trackbar);
    cv::createTrackbar("window size", "Binary", &wsize, 15, on_trackbar);
    drawWindow();

    while(true) {
        int key = cv::waitKey(0);
        if (key == 'q') {
            return 0;
        }
        if(key == 's') {
            cv::imwrite("result.png", dst);
        }
    }
}

実行結果

入力画像
f:id:kivantium:20170508230735p:plain

手動で閾値を設定してセグメンテーションした結果
f:id:kivantium:20170508230637p:plain

local intensity maximumとなるピクセルを赤く表示した結果
f:id:kivantium:20170508230819p:plain

出力は

population is: 21

でした。

もう少し難しい入力画像と処理結果
f:id:kivantium:20170508230907p:plain
f:id:kivantium:20170508231003p:plain
f:id:kivantium:20170508231012p:plain

population is: 55

これより大きい画像について実行しようとすると、画像の場所ごとに明るさが違うようで単純な閾値設定による二値化ではうまくセグメンテーションができませんでした。AdaptiveThresholdも試しましたが、OpenCVの関数ではうまくできませんでした。次回は妥協ポイントその1をなんとかする必要がありそうです。