1 /*M///////////////////////////////////////////////////////////////////////////////////////
2 //
3 // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
4 //
5 // By downloading, copying, installing or using the software you agree to this license.
6 // If you do not agree to this license, do not download, install,
7 // copy or use the software.
8 //
9 //
10 // License Agreement
11 // For Open Source Computer Vision Library
12 //
13 // Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
14 // Copyright (C) 2009, Willow Garage Inc., all rights reserved.
15 // Third party copyrights are property of their respective owners.
16 //
17 // Redistribution and use in source and binary forms, with or without modification,
18 // are permitted provided that the following conditions are met:
19 //
20 // * Redistribution's of source code must retain the above copyright notice,
21 // this list of conditions and the following disclaimer.
22 //
23 // * Redistribution's in binary form must reproduce the above copyright notice,
24 // this list of conditions and the following disclaimer in the documentation
25 // and/or other materials provided with the distribution.
26 //
27 // * The name of the copyright holders may not be used to endorse or promote products
28 // derived from this software without specific prior written permission.
29 //
30 // This software is provided by the copyright holders and contributors "as is" and
31 // any express or implied warranties, including, but not limited to, the implied
32 // warranties of merchantability and fitness for a particular purpose are disclaimed.
33 // In no event shall the Intel Corporation or contributors be liable for any direct,
34 // indirect, incidental, special, exemplary, or consequential damages
35 // (including, but not limited to, procurement of substitute goods or services;
36 // loss of use, data, or profits; or business interruption) however caused
37 // and on any theory of liability, whether in contract, strict liability,
38 // or tort (including negligence or otherwise) arising in any way out of
39 // the use of this software, even if advised of the possibility of such damage.
40 //
41 //M*/
42
43 #include "precomp.hpp"
44
45 #ifdef _MSC_VER
46 # pragma warning(disable: 4512)
47 #endif
48
49 //
50 // 2D dense optical flow algorithm from the following paper:
51 // Michael Tao, Jiamin Bai, Pushmeet Kohli, and Sylvain Paris.
52 // "SimpleFlow: A Non-iterative, Sublinear Optical Flow Algorithm"
53 // Computer Graphics Forum (Eurographics 2012)
54 // http://graphics.berkeley.edu/papers/Tao-SAN-2012-05/
55 //
56
57 namespace cv
58 {
59 namespace optflow
60 {
61
62 static const uchar MASK_TRUE_VALUE = (uchar)255;
63
dist(const Vec3b & p1,const Vec3b & p2)64 inline static int dist(const Vec3b &p1, const Vec3b &p2) {
65 int a = p1[0] - p2[0];
66 int b = p1[1] - p2[1];
67 int c = p1[2] - p2[2];
68
69 return a*a+b*b+c*c;
70 }
71
dist(const Vec2f & p1,const Vec2f & p2)72 inline static float dist(const Vec2f &p1, const Vec2f &p2) {
73 float a = p1[0] - p2[0];
74 float b = p1[1] - p2[1];
75
76 return a*a+b*b;
77 }
78
79 template<class T>
min(T t1,T t2,T t3)80 inline static T min(T t1, T t2, T t3) {
81 return (t1 <= t2 && t1 <= t3) ? t1 : min(t2, t3);
82 }
83
removeOcclusions(const Mat & flow,const Mat & flow_inv,float occ_thr,Mat & confidence)84 static void removeOcclusions(const Mat& flow,
85 const Mat& flow_inv,
86 float occ_thr,
87 Mat& confidence) {
88 const int rows = flow.rows;
89 const int cols = flow.cols;
90 if (!confidence.data) {
91 confidence = Mat::zeros(rows, cols, CV_32F);
92 }
93 for (int r = 0; r < rows; ++r) {
94 for (int c = 0; c < cols; ++c) {
95 if (dist(flow.at<Vec2f>(r, c), -flow_inv.at<Vec2f>(r, c)) > occ_thr) {
96 confidence.at<float>(r, c) = 0;
97 } else {
98 confidence.at<float>(r, c) = 1;
99 }
100 }
101 }
102 }
103
wd(Mat & d,int top_shift,int bottom_shift,int left_shift,int right_shift,double sigma)104 static void wd(Mat& d, int top_shift, int bottom_shift, int left_shift, int right_shift, double sigma) {
105 for (int dr = -top_shift, r = 0; dr <= bottom_shift; ++dr, ++r) {
106 for (int dc = -left_shift, c = 0; dc <= right_shift; ++dc, ++c) {
107 d.at<float>(r, c) = (float)-(dr*dr + dc*dc);
108 }
109 }
110 d *= 1.0 / (2.0 * sigma * sigma);
111 exp(d, d);
112 }
113
114 template<typename JointVec, typename SrcVec>
115 class CrossBilateralFilter : public ParallelLoopBody {
116 Mat &joint, &confidence, &src;
117 Mat &dst;
118 int radius;
119 bool flag;
120 Mat &spaceWeights;
121 std::vector<double> &expLut;
122
123 public:
CrossBilateralFilter(Mat & joint_,Mat & confidence_,Mat & src_,Mat & dst_,int radius_,bool flag_,Mat & spaceWeights_,std::vector<double> & expLut_)124 CrossBilateralFilter(Mat &joint_, Mat &confidence_, Mat &src_, Mat &dst_, int radius_, bool flag_, Mat &spaceWeights_, std::vector<double> &expLut_)
125 :
126 joint(joint_),
127 confidence(confidence_),
128 src(src_),
129 dst(dst_),
130 radius(radius_),
131 flag(flag_),
132 spaceWeights(spaceWeights_),
133 expLut(expLut_) {
134 CV_DbgAssert(joint.type() == traits::Type<JointVec>::value && confidence.type() == CV_32F && src.type() == dst.type() && src.type() == traits::Type<SrcVec>::value);
135 CV_DbgAssert(joint.rows == src.rows && confidence.rows == src.rows && src.rows == dst.rows + 2 * radius);
136 CV_DbgAssert(joint.cols == src.cols && confidence.cols == src.cols && src.cols == dst.cols + 2 * radius);
137 }
138
operator ()(const Range & range) const139 void operator()(const Range &range) const CV_OVERRIDE {
140 const int d = 2*radius +1;
141 for (int i = range.start; i < range.end; i++) {
142 SrcVec* dstRow = dst.ptr<SrcVec>(i);
143 for (int j = 0; j < dst.cols; j++) {
144 const JointVec& centeralPoint = joint.at<JointVec>(i+radius, j+radius);
145
146 Scalar totalSum = Scalar::all(0);
147 double weightsSum = 0;
148 for (int dr = i, r = 0; dr < i + d; ++dr, ++r) {
149 const JointVec *jointRow = joint.ptr<JointVec>(dr);
150 const SrcVec *srcRow = src.ptr<SrcVec>(dr);
151 const float *confidenceRow = confidence.ptr<float>(dr);
152 const float *spaceWeightsRow = spaceWeights.ptr<float>(r);
153 for (int dc = j, c = 0; dc < j + d; ++dc, ++c) {
154 double weight = spaceWeightsRow[c]*confidenceRow[dc];
155 for (int cn = 0; cn < JointVec::channels; cn++) {
156 weight *= expLut[std::abs(centeralPoint[cn] - jointRow[dc][cn])];
157 }
158 for (int cn = 0; cn < SrcVec::channels; cn++) {
159 totalSum[cn] += weight * srcRow[dc][cn];
160 }
161 weightsSum += weight;
162 }
163 }
164
165 SrcVec& srcSum = dstRow[j];
166 for (int cn = 0; cn < SrcVec::channels; cn++) {
167 srcSum[cn] = (flag && fabs(weightsSum) < 1e-9)
168 ? src.at<SrcVec>(i+radius, j+radius)[cn]
169 : static_cast<float>(totalSum[cn] / weightsSum);
170 }
171 }
172 }
173 }
174 };
175
crossBilateralFilter(InputArray joint_,InputArray confidence_,InputOutputArray src_,int radius,double sigmaColor,double sigmaSpace,bool flag=false)176 static void crossBilateralFilter(InputArray joint_,
177 InputArray confidence_,
178 InputOutputArray src_,
179 int radius,
180 double sigmaColor, double sigmaSpace,
181 bool flag = false) {
182 CV_Assert(!src_.empty());
183 CV_Assert(!confidence_.empty());
184 CV_Assert(!joint_.empty());
185
186 Mat src = src_.getMat();
187 Mat joint = joint_.getMat();
188 Mat confidence = confidence_.getMat();
189
190 CV_Assert(src.size() == joint.size() && confidence.size() == src.size());
191 CV_Assert(joint.depth() == CV_8U && confidence.type() == CV_32F);
192
193 if (sigmaColor <= 0)
194 sigmaColor = 1;
195 if (sigmaSpace <= 0)
196 sigmaSpace = 1;
197
198 if (radius <= 0)
199 radius = cvRound(sigmaSpace * 1.5);
200 radius = std::max(radius, 1);
201
202 if (src.data == joint.data)
203 joint = joint.clone();
204
205 int d = 2 * radius + 1;
206
207 Mat jointTemp, confidenceTemp, srcTemp;
208 copyMakeBorder(joint, jointTemp, radius, radius, radius, radius, BORDER_DEFAULT);
209 copyMakeBorder(confidence, confidenceTemp, radius, radius, radius, radius, BORDER_CONSTANT, Scalar(0));
210 copyMakeBorder(src, srcTemp, radius, radius, radius, radius, BORDER_DEFAULT);
211
212 Mat spaceWeights(d, d, CV_32F);
213 wd(spaceWeights, radius, radius, radius, radius, sigmaSpace);
214
215 double gaussColorCoeff = -0.5 / (sigmaColor * sigmaColor);
216 std::vector<double> expLut (256);
217 for (size_t i = 0; i < expLut.size(); i++) {
218 expLut[i] = std::exp(i * i * gaussColorCoeff);
219 }
220
221 Range range(0, src.rows);
222 parallel_for_(range, CrossBilateralFilter<Vec3b, Vec2f>(jointTemp, confidenceTemp, srcTemp, src, radius, flag, spaceWeights, expLut));
223 }
224
calcConfidence(const Mat & prev,const Mat & next,const Mat & flow,Mat & confidence,int max_flow)225 static void calcConfidence(const Mat& prev,
226 const Mat& next,
227 const Mat& flow,
228 Mat& confidence,
229 int max_flow) {
230 const int rows = prev.rows;
231 const int cols = prev.cols;
232 confidence = Mat::zeros(rows, cols, CV_32F);
233
234 for (int r0 = 0; r0 < rows; ++r0) {
235 for (int c0 = 0; c0 < cols; ++c0) {
236 Vec2f flow_at_point = flow.at<Vec2f>(r0, c0);
237 int u0 = cvRound(flow_at_point[0]);
238 if (r0 + u0 < 0) { u0 = -r0; }
239 if (r0 + u0 >= rows) { u0 = rows - 1 - r0; }
240 int v0 = cvRound(flow_at_point[1]);
241 if (c0 + v0 < 0) { v0 = -c0; }
242 if (c0 + v0 >= cols) { v0 = cols - 1 - c0; }
243
244 const int top_row_shift = -std::min(r0 + u0, max_flow);
245 const int bottom_row_shift = std::min(rows - 1 - (r0 + u0), max_flow);
246 const int left_col_shift = -std::min(c0 + v0, max_flow);
247 const int right_col_shift = std::min(cols - 1 - (c0 + v0), max_flow);
248
249 bool first_flow_iteration = true;
250 int sum_e = 0, min_e = 0;
251
252 for (int u = top_row_shift; u <= bottom_row_shift; ++u) {
253 for (int v = left_col_shift; v <= right_col_shift; ++v) {
254 int e = dist(prev.at<Vec3b>(r0, c0), next.at<Vec3b>(r0 + u0 + u, c0 + v0 + v));
255 if (first_flow_iteration) {
256 sum_e = e;
257 min_e = e;
258 first_flow_iteration = false;
259 } else {
260 sum_e += e;
261 min_e = std::min(min_e, e);
262 }
263 }
264 }
265 int windows_square = (bottom_row_shift - top_row_shift + 1) *
266 (right_col_shift - left_col_shift + 1);
267 confidence.at<float>(r0, c0) = (windows_square == 0) ? 0
268 : static_cast<float>(sum_e) / windows_square - min_e;
269 CV_Assert(confidence.at<float>(r0, c0) >= 0);
270 }
271 }
272 }
273
274 template<typename SrcVec, typename DstVec>
275 class CalcOpticalFlowSingleScaleSF : public ParallelLoopBody {
276 Mat &prev, &next;
277 Mat &mask;
278 Mat &dst;
279 int radius, maxFlow;
280 Mat &spaceWeights;
281 std::vector<double> &expLut;
282
283 public:
CalcOpticalFlowSingleScaleSF(Mat & prev_,Mat & next_,Mat & mask_,Mat & dst_,int radius_,int maxFlow_,Mat & spaceWeights_,std::vector<double> & expLut_)284 CalcOpticalFlowSingleScaleSF(Mat &prev_, Mat &next_, Mat &mask_, Mat &dst_, int radius_, int maxFlow_, Mat &spaceWeights_, std::vector<double> &expLut_)
285 :
286 prev(prev_),
287 next(next_),
288 mask(mask_),
289 dst(dst_),
290 radius(radius_),
291 maxFlow(maxFlow_),
292 spaceWeights(spaceWeights_),
293 expLut(expLut_){
294 CV_DbgAssert(prev.type() == next.type());
295 CV_DbgAssert(prev.rows == next.rows && prev.rows == dst.rows + 2 * radius);
296 CV_DbgAssert(prev.cols == next.cols && next.cols == dst.cols + 2 * radius);
297 }
298
operator ()(const Range & range) const299 void operator()(const Range &range) const CV_OVERRIDE {
300 int d = 2 * radius + 1;
301 Mat weights(d, d, CV_32F);
302 for (int i = range.start; i < range.end; i++) {
303 const uchar *maskRow = mask.ptr<uchar>(i);
304 const DstVec *dstRow = dst.ptr<DstVec>(i);
305 for (int j = 0; j < dst.cols; j++) {
306 if (!maskRow[j]) {
307 continue;
308 }
309
310 // TODO: do smth with this creepy staff
311 const DstVec &flowAtPoint = dstRow[j];
312 int u0 = cvRound(flowAtPoint[0]);
313 if (i + u0 < 0) {u0 = -i;}
314 if (i + u0 >= dst.rows) {u0 = dst.rows - 1 - i;}
315 int v0 = cvRound(flowAtPoint[1]);
316 if (j + v0 < 0) {v0 = -j;}
317 if (j + v0 >= dst.cols) {v0 = dst.cols - 1 - j;}
318
319 const int topRowShift = -std::min(i + u0, maxFlow);
320 const int bottomRowShift = std::min(dst.rows - 1 - (i + u0), maxFlow);
321 const int leftColShift = -std::min(j + v0, maxFlow);
322 const int rightColShift = std::min(dst.cols - 1 - (j + v0), maxFlow);
323
324 float minCost = FLT_MAX, bestU = (float) u0, bestV = (float) v0;
325
326 const SrcVec& centeralPoint = prev.at<SrcVec>(i+radius, j+radius);
327
328 int left_border = j, right_border = j + d;
329 for (int dr = i, r = 0; dr < i + d; ++dr, ++r) {
330 const SrcVec *prevRow = prev.ptr<SrcVec>(dr);
331 const float* spaceWeightsRow = spaceWeights.ptr<float>(r);
332 float *weightsRow = weights.ptr<float>(r);
333 for (int dc = left_border, c = 0; dc < right_border; ++dc, ++c) {
334 double weight = spaceWeightsRow[c];
335 for(int cn=0;cn<SrcVec::channels;cn++){
336 weight *= expLut[std::abs(centeralPoint[cn]-prevRow[dc][cn])];
337 }
338 weightsRow[c] = static_cast<float>(weight);
339 }
340 }
341
342 for (int u = topRowShift; u <= bottomRowShift; ++u) {
343 const int next_extended_top_window_row = i + u0 + u;
344 for (int v = leftColShift; v <= rightColShift; ++v) {
345 const int next_extended_left_window_col = j + v0 + v;
346
347 float cost = 0;
348 for (int r = 0; r < d; ++r) {
349 const SrcVec *prev_extended_window_row = prev.ptr<SrcVec>(i + r);
350 const SrcVec *next_extended_window_row = next.ptr<SrcVec>(next_extended_top_window_row + r);
351 const float *weight_window_row = weights.ptr<float>(r);
352 for (int c = 0; c < d; ++c) {
353 cost += weight_window_row[c] *
354 dist(prev_extended_window_row[j + c],
355 next_extended_window_row[next_extended_left_window_col + c]);
356 }
357 }
358 // cost should be divided by sum(weight_window), but because
359 // we interested only in min(cost) and sum(weight_window) is constant
360 // for every point - we remove it
361
362 if (cost < minCost) {
363 minCost = cost;
364 bestU = (float) (u + u0);
365 bestV = (float) (v + v0);
366 }
367 }
368 }
369
370 dst.at<DstVec>(i, j) = DstVec(bestU, bestV);
371 }
372 }
373 }
374 };
375
calcOpticalFlowSingleScaleSF(InputArray prev_,InputArray next_,InputArray mask_,InputOutputArray dst_,int radius,int max_flow,float sigmaSpace,float sigmaColor)376 static void calcOpticalFlowSingleScaleSF(InputArray prev_,
377 InputArray next_,
378 InputArray mask_,
379 InputOutputArray dst_,
380 int radius,
381 int max_flow,
382 float sigmaSpace,
383 float sigmaColor) {
384 Mat prev = prev_.getMat();
385 Mat next = next_.getMat();
386 Mat mask = mask_.getMat();
387 Mat dst = dst_.getMat();
388
389 Mat prevTemp, nextTemp;
390 copyMakeBorder(prev, prevTemp, radius, radius, radius, radius, BORDER_DEFAULT);
391 copyMakeBorder(next, nextTemp, radius, radius, radius, radius, BORDER_DEFAULT);
392
393 int d = 2 * radius + 1;
394
395 Mat spaceWeights(d, d, CV_32F);
396 wd(spaceWeights, radius, radius, radius, radius, sigmaSpace);
397
398 double gaussColorCoeff = -0.5 / (sigmaColor * sigmaColor);
399 std::vector<double> expLut (256);
400 for (size_t i = 0; i < expLut.size(); i++) {
401 expLut[i] = std::exp(i * i * gaussColorCoeff);
402 }
403
404 Range range(0, dst.rows);
405 parallel_for_(range, CalcOpticalFlowSingleScaleSF<Vec3b, Vec2f>(prevTemp, nextTemp, mask, dst, radius, max_flow, spaceWeights, expLut));
406 }
407
upscaleOpticalFlow(int new_rows,int new_cols,const Mat & image,const Mat & confidence,Mat & flow,int averaging_radius,float sigma_dist,float sigma_color)408 static Mat upscaleOpticalFlow(int new_rows,
409 int new_cols,
410 const Mat& image,
411 const Mat& confidence,
412 Mat& flow,
413 int averaging_radius,
414 float sigma_dist,
415 float sigma_color) {
416 crossBilateralFilter(image, confidence, flow, averaging_radius, sigma_color, sigma_dist, true);
417 Mat new_flow;
418 resize(flow, new_flow, Size(new_cols, new_rows), 0, 0, INTER_NEAREST);
419 new_flow *= 2;
420 return new_flow;
421 }
422
calcIrregularityMat(const Mat & flow,int radius)423 static Mat calcIrregularityMat(const Mat& flow, int radius) {
424 const int rows = flow.rows;
425 const int cols = flow.cols;
426 Mat irregularity = Mat::zeros(rows, cols, CV_32F);
427 for (int r = 0; r < rows; ++r) {
428 const int start_row = std::max(0, r - radius);
429 const int end_row = std::min(rows - 1, r + radius);
430 for (int c = 0; c < cols; ++c) {
431 const int start_col = std::max(0, c - radius);
432 const int end_col = std::min(cols - 1, c + radius);
433 for (int dr = start_row; dr <= end_row; ++dr) {
434 for (int dc = start_col; dc <= end_col; ++dc) {
435 const float diff = dist(flow.at<Vec2f>(r, c), flow.at<Vec2f>(dr, dc));
436 if (diff > irregularity.at<float>(r, c)) {
437 irregularity.at<float>(r, c) = diff;
438 }
439 }
440 }
441 }
442 }
443 return irregularity;
444 }
445
selectPointsToRecalcFlow(const Mat & flow,int irregularity_metric_radius,float speed_up_thr,int curr_rows,int curr_cols,const Mat & prev_speed_up,Mat & speed_up,Mat & mask)446 static void selectPointsToRecalcFlow(const Mat& flow,
447 int irregularity_metric_radius,
448 float speed_up_thr,
449 int curr_rows,
450 int curr_cols,
451 const Mat& prev_speed_up,
452 Mat& speed_up,
453 Mat& mask) {
454 const int prev_rows = flow.rows;
455 const int prev_cols = flow.cols;
456
457 Mat is_flow_regular = calcIrregularityMat(flow, irregularity_metric_radius)
458 < speed_up_thr;
459 Mat done = Mat::zeros(prev_rows, prev_cols, CV_8U);
460 speed_up = Mat::zeros(curr_rows, curr_cols, CV_8U);
461 mask = Mat::zeros(curr_rows, curr_cols, CV_8U);
462
463 for (int r = 0; r < is_flow_regular.rows; ++r) {
464 for (int c = 0; c < is_flow_regular.cols; ++c) {
465 if (!done.at<uchar>(r, c)) {
466 if (is_flow_regular.at<uchar>(r, c) &&
467 2*r + 1 < curr_rows && 2*c + 1< curr_cols) {
468
469 bool all_flow_in_region_regular = true;
470 int speed_up_at_this_point = prev_speed_up.at<uchar>(r, c);
471 int step = (1 << speed_up_at_this_point) - 1;
472 int prev_top = r;
473 int prev_bottom = std::min(r + step, prev_rows - 1);
474 int prev_left = c;
475 int prev_right = std::min(c + step, prev_cols - 1);
476
477 for (int rr = prev_top; rr <= prev_bottom; ++rr) {
478 for (int cc = prev_left; cc <= prev_right; ++cc) {
479 done.at<uchar>(rr, cc) = 1;
480 if (!is_flow_regular.at<uchar>(rr, cc)) {
481 all_flow_in_region_regular = false;
482 }
483 }
484 }
485
486 int curr_top = std::min(2 * r, curr_rows - 1);
487 int curr_bottom = std::min(2*(r + step) + 1, curr_rows - 1);
488 int curr_left = std::min(2 * c, curr_cols - 1);
489 int curr_right = std::min(2*(c + step) + 1, curr_cols - 1);
490
491 if (all_flow_in_region_regular &&
492 curr_top != curr_bottom &&
493 curr_left != curr_right) {
494 mask.at<uchar>(curr_top, curr_left) = MASK_TRUE_VALUE;
495 mask.at<uchar>(curr_bottom, curr_left) = MASK_TRUE_VALUE;
496 mask.at<uchar>(curr_top, curr_right) = MASK_TRUE_VALUE;
497 mask.at<uchar>(curr_bottom, curr_right) = MASK_TRUE_VALUE;
498 for (int rr = curr_top; rr <= curr_bottom; ++rr) {
499 for (int cc = curr_left; cc <= curr_right; ++cc) {
500 speed_up.at<uchar>(rr, cc) = (uchar)(speed_up_at_this_point + 1);
501 }
502 }
503 } else {
504 for (int rr = curr_top; rr <= curr_bottom; ++rr) {
505 for (int cc = curr_left; cc <= curr_right; ++cc) {
506 mask.at<uchar>(rr, cc) = MASK_TRUE_VALUE;
507 }
508 }
509 }
510 } else {
511 done.at<uchar>(r, c) = 1;
512 for (int dr = 0; dr <= 1; ++dr) {
513 int nr = 2*r + dr;
514 for (int dc = 0; dc <= 1; ++dc) {
515 int nc = 2*c + dc;
516 if (nr < curr_rows && nc < curr_cols) {
517 mask.at<uchar>(nr, nc) = MASK_TRUE_VALUE;
518 }
519 }
520 }
521 }
522 }
523 }
524 }
525 }
526
extrapolateValueInRect(int height,int width,float v11,float v12,float v21,float v22,int r,int c)527 static inline float extrapolateValueInRect(int height, int width,
528 float v11, float v12,
529 float v21, float v22,
530 int r, int c) {
531 if (r == 0 && c == 0) { return v11;}
532 if (r == 0 && c == width) { return v12;}
533 if (r == height && c == 0) { return v21;}
534 if (r == height && c == width) { return v22;}
535
536 CV_Assert(height > 0 && width > 0);
537 float qr = float(r) / height;
538 float pr = 1.0f - qr;
539 float qc = float(c) / width;
540 float pc = 1.0f - qc;
541
542 return v11*pr*pc + v12*pr*qc + v21*qr*pc + v22*qc*qr;
543 }
544
extrapolateFlow(Mat & flow,const Mat & speed_up)545 static void extrapolateFlow(Mat& flow,
546 const Mat& speed_up) {
547 const int rows = flow.rows;
548 const int cols = flow.cols;
549 Mat done = Mat::zeros(rows, cols, CV_8U);
550 for (int r = 0; r < rows; ++r) {
551 for (int c = 0; c < cols; ++c) {
552 if (!done.at<uchar>(r, c) && speed_up.at<uchar>(r, c) > 1) {
553 int step = (1 << speed_up.at<uchar>(r, c)) - 1;
554 int top = r;
555 int bottom = std::min(r + step, rows - 1);
556 int left = c;
557 int right = std::min(c + step, cols - 1);
558
559 int height = bottom - top;
560 int width = right - left;
561 for (int rr = top; rr <= bottom; ++rr) {
562 for (int cc = left; cc <= right; ++cc) {
563 done.at<uchar>(rr, cc) = 1;
564 Vec2f flow_at_point;
565 Vec2f top_left = flow.at<Vec2f>(top, left);
566 Vec2f top_right = flow.at<Vec2f>(top, right);
567 Vec2f bottom_left = flow.at<Vec2f>(bottom, left);
568 Vec2f bottom_right = flow.at<Vec2f>(bottom, right);
569
570 flow_at_point[0] = extrapolateValueInRect(height, width,
571 top_left[0], top_right[0],
572 bottom_left[0], bottom_right[0],
573 rr-top, cc-left);
574
575 flow_at_point[1] = extrapolateValueInRect(height, width,
576 top_left[1], top_right[1],
577 bottom_left[1], bottom_right[1],
578 rr-top, cc-left);
579 flow.at<Vec2f>(rr, cc) = flow_at_point;
580 }
581 }
582 }
583 }
584 }
585 }
586
buildPyramidWithResizeMethod(const Mat & src,std::vector<Mat> & pyramid,int layers,int interpolation_type)587 static void buildPyramidWithResizeMethod(const Mat& src,
588 std::vector<Mat>& pyramid,
589 int layers,
590 int interpolation_type) {
591 pyramid.push_back(src);
592 for (int i = 1; i <= layers; ++i) {
593 Mat prev = pyramid[i - 1];
594 if (prev.rows <= 1 || prev.cols <= 1) {
595 break;
596 }
597
598 Mat next;
599 resize(prev, next, Size((prev.cols + 1) / 2, (prev.rows + 1) / 2), 0, 0, interpolation_type);
600 pyramid.push_back(next);
601 }
602 }
603
calcOpticalFlowSF(InputArray _from,InputArray _to,OutputArray _resulted_flow,int layers,int averaging_radius,int max_flow,double sigma_dist,double sigma_color,int postprocess_window,double sigma_dist_fix,double sigma_color_fix,double occ_thr,int upscale_averaging_radius,double upscale_sigma_dist,double upscale_sigma_color,double speed_up_thr)604 CV_EXPORTS_W void calcOpticalFlowSF(InputArray _from,
605 InputArray _to,
606 OutputArray _resulted_flow,
607 int layers,
608 int averaging_radius,
609 int max_flow,
610 double sigma_dist,
611 double sigma_color,
612 int postprocess_window,
613 double sigma_dist_fix,
614 double sigma_color_fix,
615 double occ_thr,
616 int upscale_averaging_radius,
617 double upscale_sigma_dist,
618 double upscale_sigma_color,
619 double speed_up_thr)
620 {
621 Mat from = _from.getMat();
622 Mat to = _to.getMat();
623
624 std::vector<Mat> pyr_from_images;
625 std::vector<Mat> pyr_to_images;
626
627 buildPyramidWithResizeMethod(from, pyr_from_images, layers - 1, INTER_CUBIC);
628 buildPyramidWithResizeMethod(to, pyr_to_images, layers - 1, INTER_CUBIC);
629
630 CV_Assert((int)pyr_from_images.size() == layers && (int)pyr_to_images.size() == layers);
631
632 Mat curr_from, curr_to, prev_from, prev_to;
633
634 curr_from = pyr_from_images[layers - 1];
635 curr_to = pyr_to_images[layers - 1];
636
637 Mat mask = Mat::ones(curr_from.size(), CV_8U);
638 Mat mask_inv = Mat::ones(curr_from.size(), CV_8U);
639
640 Mat flow = Mat::zeros(curr_from.size(), CV_32FC2);
641 Mat flow_inv = Mat::zeros(curr_to.size(), CV_32FC2);
642
643 Mat confidence;
644 Mat confidence_inv;
645
646
647 calcOpticalFlowSingleScaleSF(curr_from,
648 curr_to,
649 mask,
650 flow,
651 averaging_radius,
652 max_flow,
653 (float)sigma_dist,
654 (float)sigma_color);
655
656 calcOpticalFlowSingleScaleSF(curr_to,
657 curr_from,
658 mask_inv,
659 flow_inv,
660 averaging_radius,
661 max_flow,
662 (float)sigma_dist,
663 (float)sigma_color);
664
665 removeOcclusions(flow,
666 flow_inv,
667 (float)occ_thr,
668 confidence);
669
670 removeOcclusions(flow_inv,
671 flow,
672 (float)occ_thr,
673 confidence_inv);
674
675 Mat speed_up = Mat::zeros(curr_from.size(), CV_8U);
676 Mat speed_up_inv = Mat::zeros(curr_from.size(), CV_8U);
677
678 for (int curr_layer = layers - 2; curr_layer >= 0; --curr_layer) {
679 curr_from = pyr_from_images[curr_layer];
680 curr_to = pyr_to_images[curr_layer];
681 prev_from = pyr_from_images[curr_layer + 1];
682 prev_to = pyr_to_images[curr_layer + 1];
683
684 const int curr_rows = curr_from.rows;
685 const int curr_cols = curr_from.cols;
686
687 Mat new_speed_up, new_speed_up_inv;
688
689 selectPointsToRecalcFlow(flow,
690 averaging_radius,
691 (float)speed_up_thr,
692 curr_rows,
693 curr_cols,
694 speed_up,
695 new_speed_up,
696 mask);
697
698 selectPointsToRecalcFlow(flow_inv,
699 averaging_radius,
700 (float)speed_up_thr,
701 curr_rows,
702 curr_cols,
703 speed_up_inv,
704 new_speed_up_inv,
705 mask_inv);
706
707 speed_up = new_speed_up;
708 speed_up_inv = new_speed_up_inv;
709
710 flow = upscaleOpticalFlow(curr_rows,
711 curr_cols,
712 prev_from,
713 confidence,
714 flow,
715 upscale_averaging_radius,
716 (float)upscale_sigma_dist,
717 (float)upscale_sigma_color);
718
719 flow_inv = upscaleOpticalFlow(curr_rows,
720 curr_cols,
721 prev_to,
722 confidence_inv,
723 flow_inv,
724 upscale_averaging_radius,
725 (float)upscale_sigma_dist,
726 (float)upscale_sigma_color);
727
728 calcConfidence(curr_from, curr_to, flow, confidence, max_flow);
729 calcOpticalFlowSingleScaleSF(curr_from,
730 curr_to,
731 mask,
732 flow,
733 averaging_radius,
734 max_flow,
735 (float)sigma_dist,
736 (float)sigma_color);
737
738 calcConfidence(curr_to, curr_from, flow_inv, confidence_inv, max_flow);
739 calcOpticalFlowSingleScaleSF(curr_to,
740 curr_from,
741 mask_inv,
742 flow_inv,
743 averaging_radius,
744 max_flow,
745 (float)sigma_dist,
746 (float)sigma_color);
747
748 extrapolateFlow(flow, speed_up);
749 extrapolateFlow(flow_inv, speed_up_inv);
750
751 //TODO: should we remove occlusions for the last stage?
752 removeOcclusions(flow, flow_inv, (float)occ_thr, confidence);
753 removeOcclusions(flow_inv, flow, (float)occ_thr, confidence_inv);
754 }
755
756 crossBilateralFilter(curr_from, confidence, flow,
757 postprocess_window, (float)sigma_color_fix, (float)sigma_dist_fix);
758
759 GaussianBlur(flow, flow, Size(3, 3), 5);
760
761 _resulted_flow.create(flow.size(), CV_32FC2);
762 Mat resulted_flow = _resulted_flow.getMat();
763 int from_to[] = {0,1 , 1,0};
764 mixChannels(&flow, 1, &resulted_flow, 1, from_to, 2);
765 }
766
calcOpticalFlowSF(InputArray from,InputArray to,OutputArray flow,int layers,int averaging_block_size,int max_flow)767 CV_EXPORTS_W void calcOpticalFlowSF(InputArray from,
768 InputArray to,
769 OutputArray flow,
770 int layers,
771 int averaging_block_size,
772 int max_flow) {
773 calcOpticalFlowSF(from, to, flow, layers, averaging_block_size, max_flow,
774 4.1, 25.5, 18, 55.0, 25.5, 0.35, 18, 55.0, 25.5, 10);
775 }
776
777 }
778 }
779
780