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

生物の実験でシャーレに生えた大腸菌のコロニーの数を数えるという作業をしたのですが、こんなものは人間の仕事ではないので自動化したいと思いました。

OpenCVのインストール

OpenCV 3.0で入った関数を使いたいので記事を書いた時点で最新版のOpenCV 3.1をビルドしてインストールします。
環境はUbuntu 15.10で、基本的にInstallation in Linux — OpenCV 2.4.12.0 documentationに従います。

sudo apt-get install build-essential
sudo apt-get install python-dev python-numpy libtbb2 libtbb-dev libjpeg-dev libpng-dev libtiff-dev libjasper-dev libdc1394-22-dev
wget https://github.com/Itseez/opencv/archive/3.1.0.zip
unzip opencv-3.1.0.zip
cd opencv-3.1.0
mkdir build
cd build
cmake -D CMAKE_BUILD_TYPE=RELEASE -D CMAKE_INSTALL_PREFIX=/usr/local ..
make -j $(nproc)
sudo make install

ドキュメントでは以上なのですが、手元の環境ではlippicvが見つからないと怒られたので、追加で

cd ..
cd 3rdparty/ippicv/unpack/ippicv_lnx/lib/intel64/libippicv.a /usr/local/lib/

を実行する必要がありました（参考: c++ - compiling code with opencv - /usr/bin/ld: cannot find -lippicv - Stack Overflow）

ラベリングによるコロニー数カウント

OpenCV 3.0からラベリングAPIが導入されたのでOpenCVを使ったラベリング · atinfinity/lab Wiki · GitHubを参考にラベリングを使ってコロニー数をカウントします。

処理の流れとしては、

adaptiveThresholdを使って二値化する
connectedComponentsを使ってラベリングし、領域数を出力する

というだけの簡単なものです。二値化のパラメータを調整できるようトラックバーもつけました。以下にソースコードを示します。

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

#include <iostream>
#include <vector>

cv::Mat src, bin, blured, dst;
// medianBlur()のパラメータkSize
int kSize = 7;
// adaptiveThresholdのパラメータblockSize
int blockSize = 7;
// adaptiveThresholdのパラメータC
int C = 2;

// コロニー数を数える関数
void countColony(void) {
    // ラベル用画像生成
    cv::Mat labelImage(src.size(), CV_32S);

    // 入力画像の平滑化
    cv::medianBlur(src, blured, kSize);
    // 平滑化画像の二値化
    cv::adaptiveThreshold(blured, bin, 255, cv::ADAPTIVE_THRESH_MEAN_C, cv::THRESH_BINARY, blockSize, C);
    // ラベリング実行
    int nLabels = cv::connectedComponents(bin, labelImage, 8);
    // ラベル数を出力
    std::cout << nLabels << std::endl;

    // ラベリング結果の描画色をランダムに選ぶ
    std::vector<cv::Vec3b> colors(nLabels);
    colors[0] = cv::Vec3b(0, 0, 0);
    for(int label = 1; label < nLabels; ++label) {
        colors[label] = cv::Vec3b((rand()&255), (rand()&255), (rand()&255));
    }

    // ラベリング結果の描画
    for(int y = 0; y < dst.rows; ++y) {
        for(int x = 0; x < dst.cols; ++x) {
            int label = labelImage.at<int>(y, x);
            cv::Vec3b &pixel = dst.at<cv::Vec3b>(y, x);
            pixel = colors[label];
        }
    }
    cv::imshow("Binary", bin);
    cv::imshow("Connected Components", dst);
}

// トラックバーに変化があった時に呼ばれる関数
void on_trackbar( int, void* ) {
    // kSizeが1より大きい奇数になることを保証
    if(kSize%2 == 0) kSize += 1;
    if(kSize < 3) kSize = 3; 
    
    // blockSizeが1より大きい奇数になることを保証
    if(blockSize%2 == 0) blockSize += 1;
    if(blockSize <= 1) blockSize = 3;

    countColony();

}

int main(int argc, const char** argv) {
    // グレースケールで画像読み込み
    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("blockSize", "Binary", &blockSize, 21, on_trackbar);
    cv::createTrackbar("C", "Binary", &C, 20, on_trackbar);
    cv::createTrackbar("kSize", "Binary", &kSize, 21, on_trackbar);

    countColony();

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

結果

大腸菌の生えたシャーレのうち、コロニー数を数えたい部分を切り抜きます。
f:id:kivantium:20160408224245p:plain:w300

g++ colonycount.cpp `pkg-config opencv --cflags --libs`
./a.out ecoli.png

のように実行すると次のようになります。
f:id:kivantium:20160408224359p:plain:w300 　 f:id:kivantium:20160408224410p:plain:w300
ほとんどのコロニーがちゃんと別の部分として認識されていることが分かります。
このプログラムでの出力は42（背景と培地が含まれるのでコロニーとしては40）、目視確認して数えたのは43だったのでそこそこの精度であることが分かります。

43くらいなら余裕で数えられるのでもっとたくさん生えている画像を与えてみました。
f:id:kivantium:20160408224714p:plain:w300
ラベリングした結果は
f:id:kivantium:20160408224740p:plain:w300
のようになりました。

多くのコロニーがつながっているため一つのものとして認識されている上に培地が複数個の領域に分裂してしまっているのでまともなカウントができていないことが分かります。connectedComponentsの仕様上どうしようもないので、この単純なアプローチで多数のコロニーをカウントするのは難しそうです。

今後の展望

今後の人生で大腸菌のコロニー数を数えることがどれくらいあるか分かりませんが、もし必要なのであれば改良版を作る必要があります。もう公開されていないiPhoneで細菌を数えるコロニーカウンタアプリのような先行研究もありましたが、手法の詳細は特に書いていないので二値化+ラベリングよりまともな方法をきちんと考える必要がありそうです。

やるとすれば