1 #include <iostream>
2 #include <fstream>
3
4 #include <opencv2/core/utility.hpp>
5 #include "opencv2/video.hpp"
6 #include "opencv2/imgcodecs.hpp"
7 #include "opencv2/highgui.hpp"
8
9 using namespace cv;
10 using namespace std;
11
isFlowCorrect(Point2f u)12 inline bool isFlowCorrect(Point2f u)
13 {
14 return !cvIsNaN(u.x) && !cvIsNaN(u.y) && fabs(u.x) < 1e9 && fabs(u.y) < 1e9;
15 }
16
computeColor(float fx,float fy)17 static Vec3b computeColor(float fx, float fy)
18 {
19 static bool first = true;
20
21 // relative lengths of color transitions:
22 // these are chosen based on perceptual similarity
23 // (e.g. one can distinguish more shades between red and yellow
24 // than between yellow and green)
25 const int RY = 15;
26 const int YG = 6;
27 const int GC = 4;
28 const int CB = 11;
29 const int BM = 13;
30 const int MR = 6;
31 const int NCOLS = RY + YG + GC + CB + BM + MR;
32 static Vec3i colorWheel[NCOLS];
33
34 if (first)
35 {
36 int k = 0;
37
38 for (int i = 0; i < RY; ++i, ++k)
39 colorWheel[k] = Vec3i(255, 255 * i / RY, 0);
40
41 for (int i = 0; i < YG; ++i, ++k)
42 colorWheel[k] = Vec3i(255 - 255 * i / YG, 255, 0);
43
44 for (int i = 0; i < GC; ++i, ++k)
45 colorWheel[k] = Vec3i(0, 255, 255 * i / GC);
46
47 for (int i = 0; i < CB; ++i, ++k)
48 colorWheel[k] = Vec3i(0, 255 - 255 * i / CB, 255);
49
50 for (int i = 0; i < BM; ++i, ++k)
51 colorWheel[k] = Vec3i(255 * i / BM, 0, 255);
52
53 for (int i = 0; i < MR; ++i, ++k)
54 colorWheel[k] = Vec3i(255, 0, 255 - 255 * i / MR);
55
56 first = false;
57 }
58
59 const float rad = sqrt(fx * fx + fy * fy);
60 const float a = atan2(-fy, -fx) / (float)CV_PI;
61
62 const float fk = (a + 1.0f) / 2.0f * (NCOLS - 1);
63 const int k0 = static_cast<int>(fk);
64 const int k1 = (k0 + 1) % NCOLS;
65 const float f = fk - k0;
66
67 Vec3b pix;
68
69 for (int b = 0; b < 3; b++)
70 {
71 const float col0 = colorWheel[k0][b] / 255.f;
72 const float col1 = colorWheel[k1][b] / 255.f;
73
74 float col = (1 - f) * col0 + f * col1;
75
76 if (rad <= 1)
77 col = 1 - rad * (1 - col); // increase saturation with radius
78 else
79 col *= .75; // out of range
80
81 pix[2 - b] = static_cast<uchar>(255.f * col);
82 }
83
84 return pix;
85 }
86
drawOpticalFlow(const Mat_<Point2f> & flow,Mat & dst,float maxmotion=-1)87 static void drawOpticalFlow(const Mat_<Point2f>& flow, Mat& dst, float maxmotion = -1)
88 {
89 dst.create(flow.size(), CV_8UC3);
90 dst.setTo(Scalar::all(0));
91
92 // determine motion range:
93 float maxrad = maxmotion;
94
95 if (maxmotion <= 0)
96 {
97 maxrad = 1;
98 for (int y = 0; y < flow.rows; ++y)
99 {
100 for (int x = 0; x < flow.cols; ++x)
101 {
102 Point2f u = flow(y, x);
103
104 if (!isFlowCorrect(u))
105 continue;
106
107 maxrad = max(maxrad, sqrt(u.x * u.x + u.y * u.y));
108 }
109 }
110 }
111
112 for (int y = 0; y < flow.rows; ++y)
113 {
114 for (int x = 0; x < flow.cols; ++x)
115 {
116 Point2f u = flow(y, x);
117
118 if (isFlowCorrect(u))
119 dst.at<Vec3b>(y, x) = computeColor(u.x / maxrad, u.y / maxrad);
120 }
121 }
122 }
123
124 // binary file format for flow data specified here:
125 // http://vision.middlebury.edu/flow/data/
writeOpticalFlowToFile(const Mat_<Point2f> & flow,const string & fileName)126 static void writeOpticalFlowToFile(const Mat_<Point2f>& flow, const string& fileName)
127 {
128 static const char FLO_TAG_STRING[] = "PIEH";
129
130 ofstream file(fileName.c_str(), ios_base::binary);
131
132 file << FLO_TAG_STRING;
133
134 file.write((const char*) &flow.cols, sizeof(int));
135 file.write((const char*) &flow.rows, sizeof(int));
136
137 for (int i = 0; i < flow.rows; ++i)
138 {
139 for (int j = 0; j < flow.cols; ++j)
140 {
141 const Point2f u = flow(i, j);
142
143 file.write((const char*) &u.x, sizeof(float));
144 file.write((const char*) &u.y, sizeof(float));
145 }
146 }
147 }
148
main(int argc,const char * argv[])149 int main(int argc, const char* argv[])
150 {
151 if (argc < 3)
152 {
153 cerr << "Usage : " << argv[0] << "<frame0> <frame1> [<output_flow>]" << endl;
154 return -1;
155 }
156
157 Mat frame0 = imread(argv[1], IMREAD_GRAYSCALE);
158 Mat frame1 = imread(argv[2], IMREAD_GRAYSCALE);
159
160 if (frame0.empty())
161 {
162 cerr << "Can't open image [" << argv[1] << "]" << endl;
163 return -1;
164 }
165 if (frame1.empty())
166 {
167 cerr << "Can't open image [" << argv[2] << "]" << endl;
168 return -1;
169 }
170
171 if (frame1.size() != frame0.size())
172 {
173 cerr << "Images should be of equal sizes" << endl;
174 return -1;
175 }
176
177 Mat_<Point2f> flow;
178 Ptr<DenseOpticalFlow> tvl1 = createOptFlow_DualTVL1();
179
180 const double start = (double)getTickCount();
181 tvl1->calc(frame0, frame1, flow);
182 const double timeSec = (getTickCount() - start) / getTickFrequency();
183 cout << "calcOpticalFlowDual_TVL1 : " << timeSec << " sec" << endl;
184
185 Mat out;
186 drawOpticalFlow(flow, out);
187
188 if (argc == 4)
189 writeOpticalFlowToFile(flow, argv[3]);
190
191 imshow("Flow", out);
192 waitKey();
193
194 return 0;
195 }
196