画像処理で大腸菌のコロニー数を数えたかった話 その2
kivantium.hateblo.jp
の続編です。
前回は先行研究を一切調べずに適当にやったので今回は少し先行研究調査をして前回よりはまともな方法を実装したいと思います。
先行研究
Google検索したところ
- 細胞画像の分析を専門にするグループがカーネギーメロン大学にあった Cell Image Analysis - CMU
- ImageJの機能を使って自動で数える方法がある https://digital.bsd.uchicago.edu/docs/cell_counting_automated_and_manual.pdf
ということが分かりました。
追記: ImageJを使うにはこの解説が分かりやすそう
satoshithermophilus.hatenablog.com
見た論文の中で一番簡単そうだったCell counting based on local intensity maxima grouping for in-situ microscopy - IEEE Conference Publicationを実装することにしました。
論文の内容と妥協ポイント
この論文では
- 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); } } }
