【OpenCV/Python】OpenCVで細胞の画像解析をやってみた

# はじめに
仕事柄、細胞の画像を撮影することが多く、細胞画像をPython版OpenCVで解析してみた。

備忘録的な意味も込めて。

今回は、培養容器接着面を培養細胞が覆った割合（細胞占有面積率、あるいはコンフルエンシー）、いわゆる「細胞のコンフル」を求めてみる。

要するに、顕微鏡画像の中の細胞の占有率を、画像解析によって数値化する。

ご意見や、もっとこうしたほうがいいなどあればコメントお願いします。

[参考にしたURL](http://labs.eecs.tottori-u.ac.jp/sd/Member/oyamada/OpenCV/html/py_tutorials/py_tutorials.html)

GitHubにこれをベースにした[プログラム](https://github.com/insilicomab/cell_confluency)を公開しています。


# 使用する画像

以下のURLの直腸平滑筋細胞の画像を使用してみました。
[コスモ・バイオ](https://www.cosmobio.co.jp/connections/p_ku_view.asp?PrimaryKeyValue=4140979&ServerKey=Primary&selPrice=1)

![SCR.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/471633/5ba088b9-d194-171b-9a45-c2e4721c25f6.jpeg)

ぱっと見、コンフル（画像中の細胞の占有率）は70～80%程度というところでしょうか。

# 必要なパッケージ
Pythonを使って画像解析するため、**OpenCV**ライブラリを読み込む。

それと、**NumPy**も読み込んでおく。

```python
# ライブラリ
import cv2
import numpy as np
```

# 画像の読み込みとグレースケール化
**imread()**関数で画像データをカラーで読み込んで、**cvtColor()**関数でグレースケール化する。

**cvtColor()**関数の第一引数は入力画像（カラー画像）。

**imread()**関数で取得したデータは BGR 形式なので、``cv2.COLOR_BGR2GRAY``を第2引数に指定する。

```python
# カラー画像の読み込み
img = cv2.imread('cell.jpg', 1)

# グレースケール化
img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
```

# 画像を二値化する
画像を二値化する。

つまり、画素値がしきい値より大きければある値（白）を割り当て、そうでなければ別の値（黒）を割り当てる。

二値化には、いろいろな方法があるらしいが、「大津の二値化」というものを使ってみた。

[画像のしきい値処理](http://labs.eecs.tottori-u.ac.jp/sd/Member/oyamada/OpenCV/html/py_tutorials/py_imgproc/py_thresholding/py_thresholding.html#thresholding)

二値化するには`threshold()`関数を使う。

`threshold()`関数の第1引数は入力画像でグレースケール画像でないといけない。
第2引数はしきい値で、画素値を識別するために使う。
第3引数は最大値でしきい値以上の値を持つ画素に対して割り当てられる値。
前述の通り、OpenCVはいくつかのしきい値処理の方法があり、第4引数にて指定する。
今回は「大津の二値化」という方法を使うため、``cv2.THRESH_BINARY+cv2.THRESH_OTSU``とする。

```python
# 大津の二値化
ret,th = cv2.threshold(img_gray,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)
```
`threshold()`関数は2つの戻り値を返す。
2つ目の戻り値``th``が、**二値化画像**となる。

# モルフォロジー変換（膨張）
画像のノイズを除去するためにモルフォロジー変換（膨張）を行う。

カーネルのサイズ（今回は5×5サイズ）に依存して物体の境界付近の全画素が黒（0）から白（1）になって消える。
カーネル内に画素値が **1** の画素が一つでも含まれれば、出力画像の注目画素の画素値を ‘1’ になる。

[モルフォロジー変換](http://labs.eecs.tottori-u.ac.jp/sd/Member/oyamada/OpenCV/html/py_tutorials/py_imgproc/py_morphological_ops/py_morphological_ops.html#morphological-ops)

```python
# カーネルの設定
kernel = np.ones((5,5),np.uint8)

# モルフォロジー変換（膨張）
th_dilation = cv2.dilate(th,kernel,iterations = 1)
```

モルフォロジー変換前（**th**）
![th.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/471633/200a7e48-312c-d146-c9bc-84bc7451a906.jpeg)




モルフォロジー変換後（**th_dilation**）
![th_dilation1.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/471633/35cd3948-1cfb-ee26-871d-fba1c15997b9.jpeg)




細胞内の黒い領域を白い領域に変換できた。

# 輪郭抽出

モルフォロジー変換によりノイズを除去した画像をもとに輪郭を抽出する。

輪郭を抽出するには`findContours()`関数を使う。

`findContours()`関数の第1引数は輪郭抽出に使う画像。
第2引数は抽出モード、第3引数は輪郭の近似方法を指定する。

`findContours()`関数の戻り値である``contours``は各輪郭の座標データがNumpyのarray形式で収められている。

``contours``を使って、`drawContours()`関数により元画像に輪郭を描画する。

全輪郭を描画する時は`drawContours()`関数の第3引数を``-1``に指定する。

[輪郭: 初めの一歩](http://labs.eecs.tottori-u.ac.jp/sd/Member/oyamada/OpenCV/html/py_tutorials/py_imgproc/py_contours/py_contours_begin/py_contours_begin.html#contours-getting-started)

```python
# 輪郭抽出
contours, hierarchy = cv2.findContours(th_dilation,
                                       cv2.RETR_LIST,
                                       cv2.CHAIN_APPROX_NONE)

# 輪郭を元画像に描画
img_contour = cv2.drawContours(img, contours, -1, (0, 255, 0), 2)
```

# 白黒面積の計算

``.size``で全体の画素数を取得する。

`countNonZero()`関数で白領域、つまり細胞領域の画素数を取得する。

全体の画素数における白領域（細胞領域）の割合を計算し、細胞の占有率を算出する。

```python
# 全体の画素数
whole_area = th_dilation.size

# 白領域の画素数
white_area = cv2.countNonZero(th_dilation)

# コンフルエンシーを表示
print('cell confluency = ' + str( int(white_area / whole_area * 100)) + ' %')
```

結果

```python
cell confluency = 75 %
```

細胞のコンフルは大体75%という結果となった。

# 画像の表示

最後に、輪郭を加えた元の画像と輪郭抽出に用いた画像を表示する。

```python
# 画像の表示
cv2.imshow('img', img)
cv2.imshow('th_dilation', th_dilation)
cv2.waitKey(0)
cv2.destroyAllWindows()
```

結果

元画像+輪郭（**img**）
![img1.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/471633/50478ad9-dfc6-59b8-262c-fd3ae53660af.jpeg)



輪郭抽出に使った画像（**th_dilation**）
![th_dilation1.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/471633/76b9925c-44be-4f4c-5b6d-701aac9f4ead.jpeg)




# 最終的なスクリプト

```python
# import packages
import cv2
import numpy as np
import matplotlib.pyplot as plt

import argparse


def run(args):
    # read the image
    img = cv2.imread(f'../data/input/{args.name}')

    # BGR to RGB
    img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
    
    # gray scale
    img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    
    # binarization using Otsu's method
    _, th = cv2.threshold(img_gray, 0, 255, cv2.THRESH_BINARY+cv2.THRESH_OTSU)
    
    # configure the kernel
    kernel = np.ones((5, 5), np.uint8)
    
    # morphological transformation(Dilation)
    th_dilation = cv2.dilate(th, kernel, iterations=1)
    
    # contour extraction
    contours, _ = cv2.findContours(th_dilation, cv2.RETR_LIST, cv2.CHAIN_APPROX_NONE)
    
    # draw the contours on the source image
    img_contour = cv2.drawContours(img.copy(), contours, -1, (0, 255, 0), 2)
    
    # total number of pixels
    whole_area = th_dilation.size
    
    # number of zero area pixels
    white_area = cv2.countNonZero(th_dilation)
    
    # calculate confluency
    confluency = int(white_area / whole_area * 100)
    
    # show confluency
    print(f'Cell Confluency: {args.name} = {confluency} %')
    
    # visualization
    fig, ax = plt.subplots(1, 3, figsize=(20, 10))
    ax[0].imshow(img)
    ax[0].set_xticks([])
    ax[0].set_yticks([])
    ax[0].text(1, 0.1, args.name, verticalalignment='top', color='red', size='x-large')
    ax[1].imshow(th_dilation, cmap='gray')
    ax[1].set_xticks([])
    ax[1].set_yticks([])
    ax[2].imshow(img_contour)
    ax[2].text(1, 0.1, f'Cell Confluency: {confluency}', verticalalignment='top', color='red', size='x-large')
    ax[2].set_xticks([])
    ax[2].set_yticks([])
    fig.savefig('../output/result.png')


if __name__ == '__main__':
    parser = argparse.ArgumentParser()
    parser.add_argument('--name', type=str)
    args = parser.parse_args()
    
    run(args)
```