移動方向に応じたオプティカルフローの可視化
細胞の移動方向に応じてグラデーションでその方向を表示したヒートマップを作成します。このページのコードでは、移動量を表すベクトルは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;
}