細胞は変形しながら移動することができます。これを細胞遊走といいます。今回は細胞の移動を可視化、定量化するためのチュートリアルです。
まず、解析の対象とする動画として、Nature誌のLeader cells regulate collective cell migration via Rac activation in the downstream signaling of integrin β1 and PI3Kという論文に含まれるものをダウンロードします。Supplementally information Videos Movie 1を使うことにします。
オプティカルフロー
オプティカルフローとは、物体の動きをベクトルで表すものです。OpenCVではいくつかのアルゴリズムが実装されています。今回は、Farnebackのアルゴリズムを実装した、calcOpticalFlowFarneback関数を使います。難しいオプティカルフローの部分はこの関数が1行で行ってくれるため、プログラムに書くのは、この関数が計算してくれたベクトルを描画するなどの周辺的なことです。オプティカルフローの計算は、2つの画像を元に行います。すなわち、画像を細かい区画に分割し、それぞれで局所的なパターンマッチのようなものを行い、そのマッチ先へのベクトルを計算します。
calcOpticalFlowFarneback(HIS_source, source, flow, 0.8, 10, 15, 3, 5, 1.1, 0);
としていますが、それぞれの引数の意味は以下の通りです。ここで、Mat型の flow はオプティカルフローの計算に一時的に必要なものと思ってください。
cv::calcOpticalFlowFarneback
void calcOpticalFlowFarneback(const Mat& prevImg, const Mat& nextImg, Mat& flow, double pyrScale, int levels, int winsize, int iterations, int polyN, double polySigma, int flags)
Gunnar Farneback のアルゴリズムを用いて，密なオプティカルフローを求めます．
prevImg – 8ビット，シングルチャンネルの1番目の入力画像．
nextImg – prevImg と同じサイズ，同じ型の2番目の入力画像．
flow – 計算されたフロー画像．サイズは prevImg と同じで，型は CV_32FC2 ．
pyrScale – それぞれの画像に対する画像ピラミッドを作るためのスケール (<1) を指定します． pyrScale=0.5 は，隣り合う各層において，各層が前の層の半分のサイズになっている古典的画像ピラミッドを意味します．
levels – 最初の画像を含む，画像ピラミッドの層の数． levels=1 は，追加の層が作成されず，元画像だけが利用されることを意味します．
winsize – 平均化窓サイズ．この値が大きくなると，画像のノイズに対するアルゴリズムの頑健さが増し，高速なモーションを検出できる場合が多くなります．しかし，ボケたモーションフィールドを生成することになります．
iterations – 画像ピラミッドの各レベルにおける，アルゴリズムの反復数．
polyN – 各ピクセルにおける多項式展開を求めるために利用される，ピクセル近傍領域のサイズ．この値が大きくなると，画像はより滑らかなサーフェイスで近似され，アルゴリズムの頑健さは増します，しかし同時に，ボケたモーションフィールドも増加します．一般的には， polyN =5 あるいは 7 となります．
polySigma – 多項式展開の基底として利用される導関数を滑らかにするための，ガウス分布の標準偏差． polyN=5 ならばpolySigma=1.1 ， polyN=7 ならば polySigma=1.5 が適当な値です．
flags –
処理フラグ．以下の値の組み合わせ：OPTFLOW_USE_INITIAL_FLOW 入力 flow をフローの初期推定値として利用します．
OPTFLOW_FARNEBACK_GAUSSIAN オプティカルフロー推定のために，同サイズのボックスフィルタの代わりに winsize*winsize サイズのガウシアンフィルタを利用します．通常，このオプションによって処理速度は低下しますが，ボックスフィルタを利用する場合よりも正確なフローを求めることができます（また，同程度の頑健性を得るためには，ガウシアン窓の winsize をより大きくしなければいけません）．
[ opencv 2.2 documentation より]
細胞遊走の動画を入力に用いると次のように各局所区画の移動方向が線分で表示されます。
#include "stdlib.h"
#include "opencv/cv.h"
#include "opencv/highgui.h"
#include "string.h"
#include "stdio.h"
#include "iostream"
#include "math.h"
using namespace std;
using namespace cv;
int main(){
VideoCapture cap("srep07656-s2.avi");
int Width=cap.get(CV_CAP_PROP_FRAME_WIDTH);
int Height=cap.get(CV_CAP_PROP_FRAME_HEIGHT);
int max_frame=cap.get(CV_CAP_PROP_FRAME_COUNT); //フレーム数
Mat source(Height,Width,CV_8UC1);
Mat HIS_source(Height,Width,CV_8UC1);
for(int frame=0; frame<max_frame; frame++){
cap>>source;
Mat disp=source.clone();
cvtColor(source, source, CV_BGR2GRAY);
if(frame>0){
vector<cv::Point2f> prev_pts;
vector<cv::Point2f> next_pts;
Size flowSize(100,100); //ベクトルの数
Point2f center = cv::Point(source.cols/2., source.rows/2.);
for(int i=0; i<flowSize.width; ++i) {
for(int j=0; j<flowSize.width; ++j) {
Point2f p(i*float(source.cols)/(flowSize.width-1),
j*float(source.rows)/(flowSize.height-1));
prev_pts.push_back((p-center)*0.95f+center);
}
}
Mat flow;
vector<float> error;
calcOpticalFlowFarneback(HIS_source, source, flow, 0.8, 10, 15, 3, 5, 1.1, 0);
// オプティカルフローの表示
std::vector<cv::Point2f>::const_iterator p = prev_pts.begin();
for(; p!=prev_pts.end(); ++p) {
const cv::Point2f& fxy = flow.at<cv::Point2f>(p->y, p->x);
cv::line(disp, *p, *p+fxy*8, cv::Scalar(0),1);
}
HIS_source=source.clone();
imshow("vector", disp);
imshow("source", HIS_source);
int c = waitKey(1);
if(c==27)return 0;
}
cout<<frame<<endl;
frame+=1;
}
waitKey(0);
return 0;
}
移動速度に応じたオプティカルフローの可視化
細胞の移動速度に応じてグラデーションでその速度を表示したヒートマップを作成します。このページのコードでは、移動量を表すベクトルはfxyというPoint型の構造体に収められています。fxyは、fxy.x または fxy.y とすることでx、yの値を取り出すことができます。したがって、ベクトルの大きさは、pow(pow(fxy.x,2)+pow(fxy.y,2),0.5) というように計算することができます。前回のコードから実際に書き換える部分は3行のみです。
細胞遊走の動画を入力に用いると次のように各局所区画の移動速度に応じて、遅いほど黒く、速いほど白く表示されます。
#include "stdlib.h"
#include "opencv/cv.h"
#include "opencv/highgui.h"
#include "string.h"
#include "stdio.h"
#include "iostream"
#include "math.h"
using namespace std;
using namespace cv;
int main(){
VideoCapture cap("/Users/YM/Desktop/srep07656-s2.avi");
int Width=cap.get(CV_CAP_PROP_FRAME_WIDTH);
int Height=cap.get(CV_CAP_PROP_FRAME_HEIGHT);
int max_frame=cap.get(CV_CAP_PROP_FRAME_COUNT); //フレーム数
Mat source(Height,Width,CV_8UC1);
Mat HIS_source(Height,Width,CV_8UC1);
for(int frame=0; frame<max_frame; frame++){
cap>>source;
Mat disp=source.clone();
cvtColor(source, source, CV_BGR2GRAY);
if(frame>0){
vector<cv::Point2f> prev_pts;
vector<cv::Point2f> next_pts;
Size flowSize(100,100);
Point2f center = cv::Point(source.cols/2., source.rows/2.);
for(int i=0; i<flowSize.width; ++i) {
for(int j=0; j<flowSize.width; ++j) {
Point2f p(i*float(source.cols)/(flowSize.width-1),
j*float(source.rows)/(flowSize.height-1));
prev_pts.push_back((p-center)*0.95f+center);
}
}
Mat flow;
vector<float> error;
calcOpticalFlowFarneback(HIS_source, source, flow, 0.8, 10, 15, 3, 5, 1.1, 0);
// オプティカルフローの表示
std::vector<cv::Point2f>::const_iterator p = prev_pts.begin();
for(; p!=prev_pts.end(); ++p) {
const cv::Point2f& fxy = flow.at<cv::Point2f>(p->y, p->x);
//cv::line(disp, *p, *p+fxy*8, cv::Scalar(0),1); //この行を削除
double length=pow(pow(fxy.x,2)+pow(fxy.y,2),0.5); //この行を追加
circle(disp, *p, 5, length*50, -1, 4); //この行を追加
}
HIS_source=source.clone();
imshow("vector", disp);
imshow("source", HIS_source);
int c = waitKey(1);
if(c==27)return 0;
}
cout<<frame<<endl;
frame+=1;
}
waitKey(0);
return 0;
}
移動方向に応じたオプティカルフローの可視化
細胞の移動方向に応じてグラデーションでその方向を表示したヒートマップを作成します。このページのコードでは、移動量を表すベクトルはfxyというPoint型の構造体に収められています。fxyは、fxy.x または fxy.y とすることでx、yの値を取り出すことができるのは移動速度に応じたヒートマップと同じです。ただし、方向をグレースケールで表現することはできないので、角度を色相とみなして表現します。ベクトルの向きは逆正接関数のatan2を使って、atan2(fxy.y, fxy.x)というように計算することができます。今回は移動する方向に応じて色相を変化させるため、途中で一度HSV色空間に変換してヒートマップを作成した後、再びBGR色空間に変換し直してから表示します。imshow関数はMatをBGR色空間のものとみなして表示するため、このように2回目の変換が必要になります。OpenCVでは0~180で色相を表すため、0,180が赤、30が黄色、60が緑、90がシアン、120が青、150がマゼンタを表すおよその色相になります。下の図はOpenCVでの色相と値の対応を示したものです。
なお、下のコードではベクトルの長さが一定以上のものに対してのみ、色をつけるようにしています。この工夫により、ノイズとなる小さな動きを無視したヒートマップが作成できます。
細胞遊走の動画を入力に用いると次のように各局所区画の移動方向に応じて異なる色がグラデーションで表示されます。
#include "stdlib.h"
#include "opencv/cv.h"
#include "opencv/highgui.h"
#include "string.h"
#include "stdio.h"
#include "iostream"
#include "math.h"
using namespace std;
using namespace cv;
#define PI 3.141592
int main(){
VideoCapture cap("srep07656-s2.avi");
int Width=cap.get(CV_CAP_PROP_FRAME_WIDTH);
int Height=cap.get(CV_CAP_PROP_FRAME_HEIGHT);
int max_frame=cap.get(CV_CAP_PROP_FRAME_COUNT); //フレーム数
Mat source(Height,Width,CV_8UC1);
Mat HIS_source(Height,Width,CV_8UC1);
for(int frame=0; frame<max_frame/2; frame++){
cap>>source;
Mat disp=source.clone();
cvtColor(source, source, CV_BGR2GRAY);
if(frame>0){
vector<cv::Point2f> prev_pts;
vector<cv::Point2f> next_pts;
Size flowSize(100,100); //ベクトルの数
Point2f center = cv::Point(source.cols/2., source.rows/2.);
for(int i=0; i<flowSize.width; ++i) {
for(int j=0; j<flowSize.width; ++j) {
Point2f p(i*float(source.cols)/(flowSize.width-1),
j*float(source.rows)/(flowSize.height-1));
prev_pts.push_back((p-center)*0.95f+center);
}
}
Mat flow;
vector<float> error;
calcOpticalFlowFarneback(HIS_source, source, flow, 0.8, 10, 15, 3, 5, 1.1, 0);
// オプティカルフローの表示
cvtColor(disp,disp,CV_BGR2HSV);
std::vector<cv::Point2f>::const_iterator p = prev_pts.begin();
for(; p!=prev_pts.end(); ++p) {
const cv::Point2f& fxy = flow.at<cv::Point2f>(p->y, p->x);
double angle = atan2(fxy.y, fxy.x);
double length = pow(pow(fxy.x,2)+pow(fxy.y,2),0.5);
if(length>1){
circle(disp, *p, 5, Scalar((angle+PI)/2/PI*180,255,255), -1, 4);
}
}
// 色相と方向の対応関係を表示
for(int a=0; a<100; a++){
double angle = (double)(a-50)/100*2*PI;
circle(disp, Point(60+50*cos(angle),60+50*sin(angle)), 5, Scalar((angle+PI)/2/PI*180,255,255), -1, 4);
}
cvtColor(disp,disp,CV_HSV2BGR);
HIS_source=source.clone();
imshow("vector", disp);
imshow("source", HIS_source);
int c = waitKey(1);
if(c==27)return 0;
}
cout<<frame<<endl;
frame+=1;
}
waitKey(0);
return 0;
}