以下の記事で画像中の円を検出する方法について考えましたが、次は特定の座標が円内にあるかないかを判定する方法について考えて見たいと思います。
目次
内外判定の対象
今回の目的は、特殊なゲル球のなかにトラップされているサイバーグリーン染色された細胞を検出することなので、染色された細胞について詳しく見ていきます。
以下の記事で閾値を手動で決定して、染色されたコロニーのみを抽出しました。
その時に取得した画像が、こちらです。
なんだかホコリみたいな画像ですね。
しかし、緑色で囲まれているものは立派なコロニーなので、一個一個を大切に扱っていきます。
判定の戦略
モーメントから重心計算
上の画像を使ってどうやってゲルの中にトラップされているかを判定するかという話ですが、まずは各コロニーの重心の座標をマッピングすることで、均一な二次元データとして扱うことができます。
import cv2
import numpy as np
img = cv2.imread('img_cells.png')
img_gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
ret, thresh = cv2.threshold(img_gray, 120, 255, cv2.THRESH_BINARY)
contours, hierarchy = cv2.findContours(thresh, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_TC89_L1)
contours_by_area = list(filter(lambda x: cv2.contourArea(x) >=0 and cv2.contourArea(x) <= 250, contours))
print(contours_by_area)
ここまでの処理で、contours_by_areaに各コロニーに対する輪郭情報が代入できているので、ループを使用して取り出していきます。
各重心の計算
変数contours_by_areaに格納されているデータのモーメントを計算し、そこから重心座標を直交座標系で表します。
ただし、今回の輪郭データには問題があります。
n>2の時は輪郭データのn次元リストからそのままモーメントを求めることができますが、n=2の時は2質点系でのモーメントの和は0になるため、モーメントから重心を計算することができません。
また、今回のコロニーの最小サイズの閾値は1ピクセルまで下げているため、1次元の輪郭情報も含まれるわけです。
そのため、三種類に場合わけして重心座標を求める必要が出てきます。
(というか、点や線には重心という概念がないですね、、)
n次元(n>2)の場合の重心計算
この場合は全ての輪郭を構成する点を考慮してモーメントを計算後、そこから重心を計算します。
cells_locations = []
for i,c in enumerate(contours_by_area):
M = cv2.moments(c)
cont = contours_by_area[i]
if len(cont) > 2:
x = M["m10"] / M["m00"]
y = M["m01"] / M["m00"]
cells_locations.append(np.array([x,y]))
2次元の場合
二次元の場合は、単純に中心座標を求めます。
elif len(cont) == 2:
x = (int(cont[0][0][0]) + int(cont[1][0][0]))/2
y = (int(cont[0][0][1]) + int(cont[1][0][1]))/2
cells_locations.append(np.array([x,y]))
ここで、輪郭を構成する座標情報は3次元データとして格納されているため、リストを参照する時にはインデックスに注意します。
1次元の場合
一次元の場合でも、輪郭データは3次元データとして保持されているので、インデックスに注意してピクセルの座標を取り出します。
elif len(cont) == 1:
x = cont[0][0][0]
y = cont[0][0][1]
cells_locations.append(np.array([x,y]))
上記の操作を行うことで、変数cells_locationsに全ての輪郭の重心を格納することができました。
hough変換のデータ構造
各ピクセルの直交座標が射影されたρ-θパラメータ空間で投票が行われ、もっともらしい円が検出されると、中心座標と半径という形で円に関するデータが返されます。
これは、内外判定には都合の良いことです。
上記で取得したコロニーの重心座標が円の中に存在するかどうかを判定するには、中学校で習った円の方程式を使います。
例えば(p,q)を中心とする半径rの円の中に重心座標(x,y)のコロニーが入っているかどうかは以下の式を満たすかどうかで判定できます。
$$(x-p)^{2}+(y-q)^{2} <= r^{2}$$
すでに、円の情報と重心座標は全て手に入っているので、あとは判定をしていくだけです。
内外判定の実装
今回使用する変数は円の情報が格納されたcirclesおよび、コロニーの重心座標が格納されたcells_locationsの二つです。
import math
from cells_contours import cells_locations
from hough_circles import circles
def validation(x,y,radius,p,q):
if ((x-p)**2 + (y-q)**2)< radius**2:
return True
if __name__ == '__main__':
matching = 0
for i in range(len(cells_locations)):
for j in range(len(circles[0])):
cell_loc_x = cells_locations[i][0]
cell_loc_y = cells_locations[i][1]
p = circles[0][j][0]
q = circles[0][j][1]
radius = circles[0][j][2]
if validation(x=cell_loc_x,y=cell_loc_y,radius=radius,p=p,q=q):
matching += 1
上記の処理では、コロニー重心座標がゲルの内側にあるかないかのみを判定して、もしコロニーの重心がゲルの内部にあった場合はカウンター変数matchingに1を代入します。
ループが2重になっていますが、そこまで重い処理ではないため画像一枚につき0.5秒ほどで終了します。