modules/video/src/simpleflow.cpp

   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 #include "simpleflow.hpp"
  45
  46 //
  47 // 2D dense optical flow algorithm from the following paper:
  48 // Michael Tao, Jiamin Bai, Pushmeet Kohli, and Sylvain Paris.
  49 // "SimpleFlow: A Non-iterative, Sublinear Optical Flow Algorithm"
  50 // Computer Graphics Forum (Eurographics 2012)
  51 // http://graphics.berkeley.edu/papers/Tao-SAN-2012-05/
  52 //
  53
  54 namespace cv
  55 {
  56
  57 static void removeOcclusions(const Mat& flow,
  58                              const Mat& flow_inv,
  59                              float occ_thr,
  60                              Mat& confidence) {
  61   const int rows = flow.rows;
  62   const int cols = flow.cols;
  63   if (!confidence.data) {
  64     confidence = Mat::zeros(rows, cols, CV_32F);
  65   }
  66   for (int r = 0; r < rows; ++r) {
  67     for (int c = 0; c < cols; ++c) {
  68       if (dist(flow.at<Vec2f>(r, c), -flow_inv.at<Vec2f>(r, c)) > occ_thr) {
  69         confidence.at<float>(r, c) = 0;
  70       } else {
  71         confidence.at<float>(r, c) = 1;
  72       }
  73     }
  74   }
  75 }
  76
  77 static void wd(Mat& d, int top_shift, int bottom_shift, int left_shift, int right_shift, float sigma) {
  78   for (int dr = -top_shift, r = 0; dr <= bottom_shift; ++dr, ++r) {
  79     for (int dc = -left_shift, c = 0; dc <= right_shift; ++dc, ++c) {
  80       d.at<float>(r, c) = (float)-(dr*dr + dc*dc);
  81     }
  82   }
  83   d *= 1.0 / (2.0 * sigma * sigma);
  84   exp(d, d);
  85 }
  86
  87 static void wc(const Mat& image, Mat& d, int r0, int c0,
  88                int top_shift, int bottom_shift, int left_shift, int right_shift, float sigma) {
  89   const Vec3b centeral_point = image.at<Vec3b>(r0, c0);
  90   int left_border = c0-left_shift, right_border = c0+right_shift;
  91   for (int dr = r0-top_shift, r = 0; dr <= r0+bottom_shift; ++dr, ++r) {
  92     const Vec3b *row = image.ptr<Vec3b>(dr);
  93     float *d_row = d.ptr<float>(r);
  94     for (int dc = left_border, c = 0; dc <= right_border; ++dc, ++c) {
  95       d_row[c] = -dist(centeral_point, row[dc]);
  96     }
  97   }
  98   d *= 1.0 / (2.0 * sigma * sigma);
  99   exp(d, d);
 100 }
 101
 102 static void crossBilateralFilter(const Mat& image,
 103                                  const Mat& edge_image,
 104                                  const Mat confidence,
 105                                  Mat& dst, int d,
 106                                  float sigma_color, float sigma_space,
 107                                  bool flag=false) {
 108   const int rows = image.rows;
 109   const int cols = image.cols;
 110   Mat image_extended, edge_image_extended, confidence_extended;
 111   copyMakeBorder(image, image_extended, d, d, d, d, BORDER_DEFAULT);
 112   copyMakeBorder(edge_image, edge_image_extended, d, d, d, d, BORDER_DEFAULT);
 113   copyMakeBorder(confidence, confidence_extended, d, d, d, d, BORDER_CONSTANT, Scalar(0));
 114   Mat weights_space(2*d+1, 2*d+1, CV_32F);
 115   wd(weights_space, d, d, d, d, sigma_space);
 116   Mat weights(2*d+1, 2*d+1, CV_32F);
 117   Mat weighted_sum(2*d+1, 2*d+1, CV_32F);
 118
 119   vector<Mat> image_extended_channels;
 120   split(image_extended, image_extended_channels);
 121
 122   for (int row = 0; row < rows; ++row) {
 123     for (int col = 0; col < cols; ++col) {
 124       wc(edge_image_extended, weights, row+d, col+d, d, d, d, d, sigma_color);
 125
 126       Range window_rows(row,row+2*d+1);
 127       Range window_cols(col,col+2*d+1);
 128
 129       multiply(weights, confidence_extended(window_rows, window_cols), weights);
 130       multiply(weights, weights_space, weights);
 131       float weights_sum = (float)sum(weights)[0];
 132
 133       for (int ch = 0; ch < 2; ++ch) {
 134         multiply(weights, image_extended_channels[ch](window_rows, window_cols), weighted_sum);
 135         float total_sum = (float)sum(weighted_sum)[0];
 136
 137         dst.at<Vec2f>(row, col)[ch] = (flag && fabs(weights_sum) < 1e-9)
 138           ? image.at<float>(row, col)
 139           : total_sum / weights_sum;
 140       }
 141     }
 142   }
 143 }
 144
 145 static void calcConfidence(const Mat& prev,
 146                            const Mat& next,
 147                            const Mat& flow,
 148                            Mat& confidence,
 149                            int max_flow) {
 150   const int rows = prev.rows;
 151   const int cols = prev.cols;
 152   confidence = Mat::zeros(rows, cols, CV_32F);
 153
 154   for (int r0 = 0; r0 < rows; ++r0) {
 155     for (int c0 = 0; c0 < cols; ++c0) {
 156       Vec2f flow_at_point = flow.at<Vec2f>(r0, c0);
 157       int u0 = cvRound(flow_at_point[0]);
 158       if (r0 + u0 < 0) { u0 = -r0; }
 159       if (r0 + u0 >= rows) { u0 = rows - 1 - r0; }
 160       int v0 = cvRound(flow_at_point[1]);
 161       if (c0 + v0 < 0) { v0 = -c0; }
 162       if (c0 + v0 >= cols) { v0 = cols - 1 - c0; }
 163
 164       const int top_row_shift = -min(r0 + u0, max_flow);
 165       const int bottom_row_shift = min(rows - 1 - (r0 + u0), max_flow);
 166       const int left_col_shift = -min(c0 + v0, max_flow);
 167       const int right_col_shift = min(cols - 1 - (c0 + v0), max_flow);
 168
 169       bool first_flow_iteration = true;
 170       float sum_e = 0, min_e = 0;
 171
 172       for (int u = top_row_shift; u <= bottom_row_shift; ++u) {
 173         for (int v = left_col_shift; v <= right_col_shift; ++v) {
 174           float e = dist(prev.at<Vec3b>(r0, c0), next.at<Vec3b>(r0 + u0 + u, c0 + v0 + v));
 175           if (first_flow_iteration) {
 176             sum_e = e;
 177             min_e = e;
 178             first_flow_iteration = false;
 179           } else {
 180             sum_e += e;
 181             min_e = std::min(min_e, e);
 182           }
 183         }
 184       }
 185       int windows_square = (bottom_row_shift - top_row_shift + 1) *
 186                            (right_col_shift - left_col_shift + 1);
 187       confidence.at<float>(r0, c0) = (windows_square == 0) ? 0
 188                                                            : sum_e / windows_square - min_e;
 189       CV_Assert(confidence.at<float>(r0, c0) >= 0);
 190     }
 191   }
 192 }
 193
 194 static void calcOpticalFlowSingleScaleSF(const Mat& prev_extended,
 195                                          const Mat& next_extended,
 196                                          const Mat& mask,
 197                                          Mat& flow,
 198                                          int averaging_radius,
 199                                          int max_flow,
 200                                          float sigma_dist,
 201                                          float sigma_color) {
 202   const int averaging_radius_2 = averaging_radius << 1;
 203   const int rows = prev_extended.rows - averaging_radius_2;
 204   const int cols = prev_extended.cols - averaging_radius_2;
 205
 206   Mat weight_window(averaging_radius_2 + 1, averaging_radius_2 + 1, CV_32F);
 207   Mat space_weight_window(averaging_radius_2 + 1, averaging_radius_2 + 1, CV_32F);
 208
 209   wd(space_weight_window, averaging_radius, averaging_radius, averaging_radius, averaging_radius, sigma_dist);
 210
 211   for (int r0 = 0; r0 < rows; ++r0) {
 212     for (int c0 = 0; c0 < cols; ++c0) {
 213       if (!mask.at<uchar>(r0, c0)) {
 214         continue;
 215       }
 216
 217       // TODO: do smth with this creepy staff
 218       Vec2f flow_at_point = flow.at<Vec2f>(r0, c0);
 219       int u0 = cvRound(flow_at_point[0]);
 220       if (r0 + u0 < 0) { u0 = -r0; }
 221       if (r0 + u0 >= rows) { u0 = rows - 1 - r0; }
 222       int v0 = cvRound(flow_at_point[1]);
 223       if (c0 + v0 < 0) { v0 = -c0; }
 224       if (c0 + v0 >= cols) { v0 = cols - 1 - c0; }
 225
 226       const int top_row_shift = -min(r0 + u0, max_flow);
 227       const int bottom_row_shift = min(rows - 1 - (r0 + u0), max_flow);
 228       const int left_col_shift = -min(c0 + v0, max_flow);
 229       const int right_col_shift = min(cols - 1 - (c0 + v0), max_flow);
 230
 231       float min_cost = FLT_MAX, best_u = (float)u0, best_v = (float)v0;
 232
 233       wc(prev_extended, weight_window, r0 + averaging_radius, c0 + averaging_radius,
 234          averaging_radius, averaging_radius, averaging_radius, averaging_radius, sigma_color);
 235       multiply(weight_window, space_weight_window, weight_window);
 236
 237       const int prev_extended_top_window_row = r0;
 238       const int prev_extended_left_window_col = c0;
 239
 240       for (int u = top_row_shift; u <= bottom_row_shift; ++u) {
 241         const int next_extended_top_window_row = r0 + u0 + u;
 242         for (int v = left_col_shift; v <= right_col_shift; ++v) {
 243           const int next_extended_left_window_col = c0 + v0 + v;
 244
 245           float cost = 0;
 246           for (int r = 0; r <= averaging_radius_2; ++r) {
 247             const Vec3b *prev_extended_window_row = prev_extended.ptr<Vec3b>(prev_extended_top_window_row + r);
 248             const Vec3b *next_extended_window_row = next_extended.ptr<Vec3b>(next_extended_top_window_row + r);
 249             const float* weight_window_row = weight_window.ptr<float>(r);
 250             for (int c = 0; c <= averaging_radius_2; ++c) {
 251               cost += weight_window_row[c] *
 252                            dist(prev_extended_window_row[prev_extended_left_window_col + c],
 253                                 next_extended_window_row[next_extended_left_window_col + c]);
 254             }
 255           }
 256           // cost should be divided by sum(weight_window), but because
 257           // we interested only in min(cost) and sum(weight_window) is constant
 258           // for every point - we remove it
 259
 260           if (cost < min_cost) {
 261             min_cost = cost;
 262             best_u = (float)(u + u0);
 263             best_v = (float)(v + v0);
 264           }
 265         }
 266       }
 267       flow.at<Vec2f>(r0, c0) = Vec2f(best_u, best_v);
 268     }
 269   }
 270 }
 271
 272 static Mat upscaleOpticalFlow(int new_rows,
 273                                int new_cols,
 274                                const Mat& image,
 275                                const Mat& confidence,
 276                                Mat& flow,
 277                                int averaging_radius,
 278                                float sigma_dist,
 279                                float sigma_color) {
 280   crossBilateralFilter(flow, image, confidence, flow, averaging_radius, sigma_color, sigma_dist, true);
 281   Mat new_flow;
 282   resize(flow, new_flow, Size(new_cols, new_rows), 0, 0, INTER_NEAREST);
 283   new_flow *= 2;
 284   return new_flow;
 285 }
 286
 287 static Mat calcIrregularityMat(const Mat& flow, int radius) {
 288   const int rows = flow.rows;
 289   const int cols = flow.cols;
 290   Mat irregularity(rows, cols, CV_32F);
 291   for (int r = 0; r < rows; ++r) {
 292     const int start_row = max(0, r - radius);
 293     const int end_row = min(rows - 1, r + radius);
 294     for (int c = 0; c < cols; ++c) {
 295       const int start_col = max(0, c - radius);
 296       const int end_col = min(cols - 1, c + radius);
 297       for (int dr = start_row; dr <= end_row; ++dr) {
 298         for (int dc = start_col; dc <= end_col; ++dc) {
 299           const float diff = dist(flow.at<Vec2f>(r, c), flow.at<Vec2f>(dr, dc));
 300           if (diff > irregularity.at<float>(r, c)) {
 301             irregularity.at<float>(r, c) = diff;
 302           }
 303         }
 304       }
 305     }
 306   }
 307   return irregularity;
 308 }
 309
 310 static void selectPointsToRecalcFlow(const Mat& flow,
 311                                      int irregularity_metric_radius,
 312                                      float speed_up_thr,
 313                                      int curr_rows,
 314                                      int curr_cols,
 315                                      const Mat& prev_speed_up,
 316                                      Mat& speed_up,
 317                                      Mat& mask) {
 318   const int prev_rows = flow.rows;
 319   const int prev_cols = flow.cols;
 320
 321   Mat is_flow_regular = calcIrregularityMat(flow, irregularity_metric_radius)
 322                               < speed_up_thr;
 323   Mat done = Mat::zeros(prev_rows, prev_cols, CV_8U);
 324   speed_up = Mat::zeros(curr_rows, curr_cols, CV_8U);
 325   mask = Mat::zeros(curr_rows, curr_cols, CV_8U);
 326
 327   for (int r = 0; r < is_flow_regular.rows; ++r) {
 328     for (int c = 0; c < is_flow_regular.cols; ++c) {
 329       if (!done.at<uchar>(r, c)) {
 330         if (is_flow_regular.at<uchar>(r, c) &&
 331             2*r + 1 < curr_rows && 2*c + 1< curr_cols) {
 332
 333           bool all_flow_in_region_regular = true;
 334           int speed_up_at_this_point = prev_speed_up.at<uchar>(r, c);
 335           int step = (1 << speed_up_at_this_point) - 1;
 336           int prev_top = r;
 337           int prev_bottom = std::min(r + step, prev_rows - 1);
 338           int prev_left = c;
 339           int prev_right = std::min(c + step, prev_cols - 1);
 340
 341           for (int rr = prev_top; rr <= prev_bottom; ++rr) {
 342             for (int cc = prev_left; cc <= prev_right; ++cc) {
 343               done.at<uchar>(rr, cc) = 1;
 344               if (!is_flow_regular.at<uchar>(rr, cc)) {
 345                 all_flow_in_region_regular = false;
 346               }
 347             }
 348           }
 349
 350           int curr_top = std::min(2 * r, curr_rows - 1);
 351           int curr_bottom = std::min(2*(r + step) + 1, curr_rows - 1);
 352           int curr_left = std::min(2 * c, curr_cols - 1);
 353           int curr_right = std::min(2*(c + step) + 1, curr_cols - 1);
 354
 355           if (all_flow_in_region_regular &&
 356               curr_top != curr_bottom &&
 357               curr_left != curr_right) {
 358             mask.at<uchar>(curr_top, curr_left) = MASK_TRUE_VALUE;
 359             mask.at<uchar>(curr_bottom, curr_left) = MASK_TRUE_VALUE;
 360             mask.at<uchar>(curr_top, curr_right) = MASK_TRUE_VALUE;
 361             mask.at<uchar>(curr_bottom, curr_right) = MASK_TRUE_VALUE;
 362             for (int rr = curr_top; rr <= curr_bottom; ++rr) {
 363               for (int cc = curr_left; cc <= curr_right; ++cc) {
 364                 speed_up.at<uchar>(rr, cc) = (uchar)(speed_up_at_this_point + 1);
 365               }
 366             }
 367           } else {
 368             for (int rr = curr_top; rr <= curr_bottom; ++rr) {
 369               for (int cc = curr_left; cc <= curr_right; ++cc) {
 370                 mask.at<uchar>(rr, cc) = MASK_TRUE_VALUE;
 371               }
 372             }
 373           }
 374         } else {
 375           done.at<uchar>(r, c) = 1;
 376           for (int dr = 0; dr <= 1; ++dr) {
 377             int nr = 2*r + dr;
 378             for (int dc = 0; dc <= 1; ++dc) {
 379               int nc = 2*c + dc;
 380               if (nr < curr_rows && nc < curr_cols) {
 381                 mask.at<uchar>(nr, nc) = MASK_TRUE_VALUE;
 382               }
 383             }
 384           }
 385         }
 386       }
 387     }
 388   }
 389 }
 390
 391 static inline float extrapolateValueInRect(int height, int width,
 392                                             float v11, float v12,
 393                                             float v21, float v22,
 394                                             int r, int c) {
 395   if (r == 0 && c == 0) { return v11;}
 396   if (r == 0 && c == width) { return v12;}
 397   if (r == height && c == 0) { return v21;}
 398   if (r == height && c == width) { return v22;}
 399
 400   float qr = float(r) / height;
 401   float pr = 1.0f - qr;
 402   float qc = float(c) / width;
 403   float pc = 1.0f - qc;
 404
 405   return v11*pr*pc + v12*pr*qc + v21*qr*pc + v22*qc*qr;
 406 }
 407
 408 static void extrapolateFlow(Mat& flow,
 409                             const Mat& speed_up) {
 410   const int rows = flow.rows;
 411   const int cols = flow.cols;
 412   Mat done(rows, cols, CV_8U);
 413   for (int r = 0; r < rows; ++r) {
 414     for (int c = 0; c < cols; ++c) {
 415       if (!done.at<uchar>(r, c) && speed_up.at<uchar>(r, c) > 1) {
 416         int step = (1 << speed_up.at<uchar>(r, c)) - 1;
 417         int top = r;
 418         int bottom = std::min(r + step, rows - 1);
 419         int left = c;
 420         int right = std::min(c + step, cols - 1);
 421
 422         int height = bottom - top;
 423         int width = right - left;
 424         for (int rr = top; rr <= bottom; ++rr) {
 425           for (int cc = left; cc <= right; ++cc) {
 426             done.at<uchar>(rr, cc) = 1;
 427             Vec2f flow_at_point;
 428             Vec2f top_left = flow.at<Vec2f>(top, left);
 429             Vec2f top_right = flow.at<Vec2f>(top, right);
 430             Vec2f bottom_left = flow.at<Vec2f>(bottom, left);
 431             Vec2f bottom_right = flow.at<Vec2f>(bottom, right);
 432
 433             flow_at_point[0] = extrapolateValueInRect(height, width,
 434                                                       top_left[0], top_right[0],
 435                                                       bottom_left[0], bottom_right[0],
 436                                                       rr-top, cc-left);
 437
 438             flow_at_point[1] = extrapolateValueInRect(height, width,
 439                                                       top_left[1], top_right[1],
 440                                                       bottom_left[1], bottom_right[1],
 441                                                       rr-top, cc-left);
 442             flow.at<Vec2f>(rr, cc) = flow_at_point;
 443           }
 444         }
 445       }
 446     }
 447   }
 448 }
 449
 450 static void buildPyramidWithResizeMethod(Mat& src,
 451                                   vector<Mat>& pyramid,
 452                                   int layers,
 453                                   int interpolation_type) {
 454   pyramid.push_back(src);
 455   for (int i = 1; i <= layers; ++i) {
 456     Mat prev = pyramid[i - 1];
 457     if (prev.rows <= 1 || prev.cols <= 1) {
 458       break;
 459     }
 460
 461     Mat next;
 462     resize(prev, next, Size((prev.cols + 1) / 2, (prev.rows + 1) / 2), 0, 0, interpolation_type);
 463     pyramid.push_back(next);
 464   }
 465 }
 466
 467 CV_EXPORTS_W void calcOpticalFlowSF(Mat& from,
 468                                     Mat& to,
 469                                     Mat& resulted_flow,
 470                                     int layers,
 471                                     int averaging_radius,
 472                                     int max_flow,
 473                                     double sigma_dist,
 474                                     double sigma_color,
 475                                     int postprocess_window,
 476                                     double sigma_dist_fix,
 477                                     double sigma_color_fix,
 478                                     double occ_thr,
 479                                     int upscale_averaging_radius,
 480                                     double upscale_sigma_dist,
 481                                     double upscale_sigma_color,
 482                                     double speed_up_thr) {
 483   vector<Mat> pyr_from_images;
 484   vector<Mat> pyr_to_images;
 485
 486   buildPyramidWithResizeMethod(from, pyr_from_images, layers - 1, INTER_CUBIC);
 487   buildPyramidWithResizeMethod(to, pyr_to_images, layers - 1, INTER_CUBIC);
 488
 489   CV_Assert((int)pyr_from_images.size() == layers && (int)pyr_to_images.size() == layers);
 490
 491   Mat curr_from, curr_to, prev_from, prev_to;
 492   Mat curr_from_extended, curr_to_extended;
 493
 494   curr_from = pyr_from_images[layers - 1];
 495   curr_to = pyr_to_images[layers - 1];
 496
 497   copyMakeBorder(curr_from, curr_from_extended,
 498                  averaging_radius, averaging_radius, averaging_radius, averaging_radius,
 499                  BORDER_DEFAULT);
 500   copyMakeBorder(curr_to, curr_to_extended,
 501                  averaging_radius, averaging_radius, averaging_radius, averaging_radius,
 502                  BORDER_DEFAULT);
 503
 504   Mat mask = Mat::ones(curr_from.size(), CV_8U);
 505   Mat mask_inv = Mat::ones(curr_from.size(), CV_8U);
 506
 507   Mat flow(curr_from.size(), CV_32FC2);
 508   Mat flow_inv(curr_to.size(), CV_32FC2);
 509
 510   Mat confidence;
 511   Mat confidence_inv;
 512
 513
 514   calcOpticalFlowSingleScaleSF(curr_from_extended,
 515                                curr_to_extended,
 516                                mask,
 517                                flow,
 518                                averaging_radius,
 519                                max_flow,
 520                                (float)sigma_dist,
 521                                (float)sigma_color);
 522
 523   calcOpticalFlowSingleScaleSF(curr_to_extended,
 524                                curr_from_extended,
 525                                mask_inv,
 526                                flow_inv,
 527                                averaging_radius,
 528                                max_flow,
 529                                (float)sigma_dist,
 530                                (float)sigma_color);
 531
 532   removeOcclusions(flow,
 533                    flow_inv,
 534                    (float)occ_thr,
 535                    confidence);
 536
 537   removeOcclusions(flow_inv,
 538                    flow,
 539                    (float)occ_thr,
 540                    confidence_inv);
 541
 542   Mat speed_up = Mat::zeros(curr_from.size(), CV_8U);
 543   Mat speed_up_inv = Mat::zeros(curr_from.size(), CV_8U);
 544
 545   for (int curr_layer = layers - 2; curr_layer >= 0; --curr_layer) {
 546     curr_from = pyr_from_images[curr_layer];
 547     curr_to = pyr_to_images[curr_layer];
 548     prev_from = pyr_from_images[curr_layer + 1];
 549     prev_to = pyr_to_images[curr_layer + 1];
 550
 551     copyMakeBorder(curr_from, curr_from_extended,
 552                    averaging_radius, averaging_radius, averaging_radius, averaging_radius,
 553                    BORDER_DEFAULT);
 554     copyMakeBorder(curr_to, curr_to_extended,
 555                    averaging_radius, averaging_radius, averaging_radius, averaging_radius,
 556                    BORDER_DEFAULT);
 557
 558     const int curr_rows = curr_from.rows;
 559     const int curr_cols = curr_from.cols;
 560
 561     Mat new_speed_up, new_speed_up_inv;
 562
 563     selectPointsToRecalcFlow(flow,
 564                              averaging_radius,
 565                              (float)speed_up_thr,
 566                              curr_rows,
 567                              curr_cols,
 568                              speed_up,
 569                              new_speed_up,
 570                              mask);
 571
 572     selectPointsToRecalcFlow(flow_inv,
 573                              averaging_radius,
 574                              (float)speed_up_thr,
 575                              curr_rows,
 576                              curr_cols,
 577                              speed_up_inv,
 578                              new_speed_up_inv,
 579                              mask_inv);
 580
 581     speed_up = new_speed_up;
 582     speed_up_inv = new_speed_up_inv;
 583
 584     flow = upscaleOpticalFlow(curr_rows,
 585                               curr_cols,
 586                               prev_from,
 587                               confidence,
 588                               flow,
 589                               upscale_averaging_radius,
 590                               (float)upscale_sigma_dist,
 591                               (float)upscale_sigma_color);
 592
 593     flow_inv = upscaleOpticalFlow(curr_rows,
 594                                   curr_cols,
 595                                   prev_to,
 596                                   confidence_inv,
 597                                   flow_inv,
 598                                   upscale_averaging_radius,
 599                                   (float)upscale_sigma_dist,
 600                                   (float)upscale_sigma_color);
 601
 602     calcConfidence(curr_from, curr_to, flow, confidence, max_flow);
 603     calcOpticalFlowSingleScaleSF(curr_from_extended,
 604                                  curr_to_extended,
 605                                  mask,
 606                                  flow,
 607                                  averaging_radius,
 608                                  max_flow,
 609                                  (float)sigma_dist,
 610                                  (float)sigma_color);
 611
 612     calcConfidence(curr_to, curr_from, flow_inv, confidence_inv, max_flow);
 613     calcOpticalFlowSingleScaleSF(curr_to_extended,
 614                                  curr_from_extended,
 615                                  mask_inv,
 616                                  flow_inv,
 617                                  averaging_radius,
 618                                  max_flow,
 619                                  (float)sigma_dist,
 620                                  (float)sigma_color);
 621
 622     extrapolateFlow(flow, speed_up);
 623     extrapolateFlow(flow_inv, speed_up_inv);
 624
 625     //TODO: should we remove occlusions for the last stage?
 626     removeOcclusions(flow, flow_inv, (float)occ_thr, confidence);
 627     removeOcclusions(flow_inv, flow, (float)occ_thr, confidence_inv);
 628   }
 629
 630   crossBilateralFilter(flow, curr_from, confidence, flow,
 631                        postprocess_window, (float)sigma_color_fix, (float)sigma_dist_fix);
 632
 633   GaussianBlur(flow, flow, Size(3, 3), 5);
 634
 635   resulted_flow = Mat(flow.size(), CV_32FC2);
 636   int from_to[] = {0,1 , 1,0};
 637   mixChannels(&flow, 1, &resulted_flow, 1, from_to, 2);
 638 }
 639
 640 CV_EXPORTS_W void calcOpticalFlowSF(Mat& from,
 641                                     Mat& to,
 642                                     Mat& flow,
 643                                     int layers,
 644                                     int averaging_block_size,
 645                                     int max_flow) {
 646   calcOpticalFlowSF(from, to, flow, layers, averaging_block_size, max_flow,
 647                     4.1, 25.5, 18, 55.0, 25.5, 0.35, 18, 55.0, 25.5, 10);
 648 }
 649
 650 }
 651