Updates for SimpleFlow algorithm
authorYury Zemlyanskiy <yuri.zemlyanskiy@gmail.com>
Mon, 27 Aug 2012 10:40:57 +0000 (14:40 +0400)
committerVadim Pisarevsky <vadim.pisarevsky@itseez.com>
Fri, 31 Aug 2012 10:39:59 +0000 (14:39 +0400)
+ New format for flow data - CV_32C2
+ Memory optimization
+ Cross Bilateral Filter optimization
+ Minor optimizations
+ Sample for calcOpticalFlowSF improved

modules/video/include/opencv2/video/tracking.hpp
modules/video/src/simpleflow.cpp
modules/video/src/simpleflow.hpp
modules/video/test/test_simpleflow.cpp
samples/cpp/simpleflow_demo.cpp

index 85c1881..6800c63 100644 (file)
@@ -329,9 +329,8 @@ CV_EXPORTS_W Mat estimateRigidTransform( InputArray src, InputArray dst,
   
 //! computes dense optical flow using Simple Flow algorithm
 CV_EXPORTS_W void calcOpticalFlowSF(Mat& from, 
-                                    Mat& to, 
-                                    Mat& flowX,
-                                    Mat& flowY,
+                                    Mat& to,
+                                    Mat& flow,
                                     int layers,
                                     int averaging_block_size, 
                                     int max_flow,
index 1fda361..0cd320d 100644 (file)
 namespace cv
 {
 
-WeightedCrossBilateralFilter::WeightedCrossBilateralFilter(
-  const Mat& _image,
-  int _windowSize,
-  double _sigmaDist,
-  double _sigmaColor)
-  : image(_image), 
-    windowSize(_windowSize), 
-    sigmaDist(_sigmaDist), 
-    sigmaColor(_sigmaColor) {
-    
-  expDist.resize(2*windowSize*windowSize+1);
-  const double sigmaDistSqr = 2 * sigmaDist * sigmaDist;
-  for (int i = 0; i <= 2*windowSize*windowSize; ++i) {
-    expDist[i] = exp(-i/sigmaDistSqr);
-  }
-
-  const double sigmaColorSqr = 2 * sigmaColor * sigmaColor;
-  wc.resize(image.rows);
-  for (int row = 0; row < image.rows; ++row) {
-    wc[row].resize(image.cols);
-    for (int col = 0; col < image.cols; ++col) {
-      int beginRow = max(0, row - windowSize);
-      int beginCol = max(0, col - windowSize);
-      int endRow = min(image.rows - 1, row + windowSize);
-      int endCol = min(image.cols - 1, col + windowSize);
-      wc[row][col] = build<double>(endRow - beginRow + 1, endCol - beginCol + 1);
-
-      for (int r = beginRow; r <= endRow; ++r) {
-        for (int c = beginCol; c <= endCol; ++c) {
-          wc[row][col][r - beginRow][c - beginCol] = 
-            exp(-dist(image.at<Vec3b>(row, col), 
-                      image.at<Vec3b>(r, c)) 
-                / sigmaColorSqr);
-        }
-      }
-    }
-  }
-}
-
-Mat WeightedCrossBilateralFilter::apply(Mat& matrix, Mat& weights) {
-  int rows = matrix.rows;
-  int cols = matrix.cols;
-
-  Mat result = Mat::zeros(rows, cols, CV_64F);
-  for (int row = 0; row < rows; ++row) {
-    for(int col = 0; col < cols; ++col) {
-      result.at<double>(row, col) = 
-        convolution(matrix, row, col, weights);
-    }
-  }
-  return result;
-}
-
-double WeightedCrossBilateralFilter::convolution(Mat& matrix, 
-                                                 int row, int col, 
-                                                 Mat& weights) {
-  double result = 0, weightsSum = 0;
-  int beginRow = max(0, row - windowSize);
-  int beginCol = max(0, col - windowSize);
-  int endRow = min(matrix.rows - 1, row + windowSize);
-  int endCol = min(matrix.cols - 1, col + windowSize);
-  for (int r = beginRow; r <= endRow; ++r) {
-    double* ptr = matrix.ptr<double>(r);
-    for (int c = beginCol; c <= endCol; ++c) {
-      const double w = expDist[dist(row, col, r, c)] *
-                       wc[row][col][r - beginRow][c - beginCol] *
-                       weights.at<double>(r, c);
-      result += ptr[c] * w;
-      weightsSum += w;
-    }
-  }
-  return result / weightsSum;
-}
-
-static void removeOcclusions(const Flow& flow, 
-                             const Flow& flow_inv,
-                             double occ_thr,
+static void removeOcclusions(const Mat& flow, 
+                             const Mat& flow_inv,
+                             float occ_thr,
                              Mat& confidence) {
-  const int rows = flow.u.rows;
-  const int cols = flow.v.cols;
-  int occlusions = 0;
+  const int rows = flow.rows;
+  const int cols = flow.cols;
   for (int r = 0; r < rows; ++r) {
     for (int c = 0; c < cols; ++c) {
-      if (dist(flow.u.at<double>(r, c), flow.v.at<double>(r, c),
-               -flow_inv.u.at<double>(r, c), -flow_inv.v.at<double>(r, c)) > occ_thr) {
-        confidence.at<double>(r, c) = 0;
-        occlusions++;
+      if (dist(flow.at<Vec2f>(r, c), -flow_inv.at<Vec2f>(r, c)) > occ_thr) {
+        confidence.at<float>(r, c) = 0;
+      } else {
+        confidence.at<float>(r, c) = 1;
       }
     }
   }
 }
 
-static Mat wd(int top_shift, int bottom_shift, int left_shift, int right_shift, double sigma) {
-  const double factor = 1.0 / (2.0 * sigma * sigma);
-  Mat d = Mat(top_shift + bottom_shift + 1, right_shift + left_shift + 1, CV_64F);
+static void wd(Mat& d, int top_shift, int bottom_shift, int left_shift, int right_shift, float sigma) {
+  const float factor = 1.0 / (2.0 * sigma * sigma);
   for (int dr = -top_shift, r = 0; dr <= bottom_shift; ++dr, ++r) {
     for (int dc = -left_shift, c = 0; dc <= right_shift; ++dc, ++c) {
-      d.at<double>(r, c) = -(dr*dr + dc*dc) * factor;
+      d.at<float>(r, c) = -(dr*dr + dc*dc) * factor;
     }
   }
-  Mat ed;
-  exp(d, ed);
-  return ed;
+  exp(d, d);
 }
 
-static Mat wc(const Mat& image, int r0, int c0, int top_shift, int bottom_shift, int left_shift, int right_shift, double sigma) {
-  const double factor = 1.0 / (2.0 * sigma * sigma);
-  Mat d = Mat(top_shift + bottom_shift + 1, right_shift + left_shift + 1, CV_64F);
+static void wc(const Mat& image, Mat& d, int r0, int c0, 
+               int top_shift, int bottom_shift, int left_shift, int right_shift, float sigma) {
+  const float factor = 1.0 / (2.0 * sigma * sigma);
+  const Vec3b centeral_point = image.at<Vec3b>(r0, c0);
   for (int dr = r0-top_shift, r = 0; dr <= r0+bottom_shift; ++dr, ++r) {
+    const Vec3b *row = image.ptr<Vec3b>(dr); 
+    float *d_row = d.ptr<float>(r);
     for (int dc = c0-left_shift, c = 0; dc <= c0+right_shift; ++dc, ++c) {
-      d.at<double>(r, c) = -dist(image.at<Vec3b>(r0, c0), image.at<Vec3b>(dr, dc)) * factor;
+      d_row[c] = -dist(centeral_point, row[dc]) * factor;
     }
   }
-  Mat ed;
-  exp(d, ed);
-  return ed;
+  exp(d, d);
 }
 
-inline static void dist(const Mat& m1, const Mat& m2, Mat& result) {
+static void dist(const Mat& m1, const Mat& m2, Mat& result) {
   const int rows = m1.rows;
   const int cols = m1.cols;
   for (int r = 0; r < rows; ++r) {
     const Vec3b *m1_row = m1.ptr<Vec3b>(r);
     const Vec3b *m2_row = m2.ptr<Vec3b>(r);
-    double* row = result.ptr<double>(r);
+    float* row = result.ptr<float>(r);
     for (int c = 0; c < cols; ++c) {
       row[c] = dist(m1_row[c], m2_row[c]);
     }
   }
 }
 
+static void crossBilateralFilter(const Mat& image, const Mat& edge_image, const Mat confidence, Mat& dst, int d, float sigma_color, float sigma_space, bool flag=false) {
+  const int rows = image.rows;
+  const int cols = image.cols;
+  Mat image_extended, edge_image_extended, confidence_extended;
+  copyMakeBorder(image, image_extended, d, d, d, d, BORDER_DEFAULT);
+  copyMakeBorder(edge_image, edge_image_extended, d, d, d, d, BORDER_DEFAULT);
+  copyMakeBorder(confidence, confidence_extended, d, d, d, d, BORDER_CONSTANT, Scalar(0));
+  Mat weights_space(2*d+1, 2*d+1, CV_32F);
+  wd(weights_space, d, d, d, d, sigma_space);
+  Mat weights(2*d+1, 2*d+1, CV_32F);
+  Mat weighted_sum(2*d+1, 2*d+1, CV_32F);
+  
+
+  vector<Mat> image_extended_channels;
+  split(image_extended, image_extended_channels);
+
+  for (int row = 0; row < rows; ++row) {
+    for (int col = 0; col < cols; ++col) {
+      wc(edge_image_extended, weights, row+d, col+d, d, d, d, d, sigma_color);
+
+      Range window_rows(row,row+2*d+1);
+      Range window_cols(col,col+2*d+1);
+
+      multiply(weights, confidence_extended(window_rows, window_cols), weights);
+      multiply(weights, weights_space, weights);
+      float weights_sum = sum(weights)[0];
+
+      for (int ch = 0; ch < 2; ++ch) {
+        multiply(weights, image_extended_channels[ch](window_rows, window_cols), weighted_sum);
+        float total_sum = sum(weighted_sum)[0];
+
+        dst.at<Vec2f>(row, col)[ch] = (flag && fabs(weights_sum) < 1e-9) 
+          ? image.at<float>(row, col) 
+          : total_sum / weights_sum;
+      }
+    }
+  }
+}
+
 static void calcOpticalFlowSingleScaleSF(const Mat& prev, 
                                          const Mat& next,
                                          const Mat& mask,
-                                         Flow& flow,
+                                         Mat& flow,
                                          Mat& confidence,
                                          int averaging_radius, 
                                          int max_flow,
-                                         double sigma_dist,
-                                         double sigma_color) {
+                                         float sigma_dist,
+                                         float sigma_color) {
   const int rows = prev.rows;
   const int cols = prev.cols;
-  confidence = Mat::zeros(rows, cols, CV_64F);
+  confidence = Mat::zeros(rows, cols, CV_32F);
+  
+  Mat diff_storage(averaging_radius*2 + 1, averaging_radius*2 + 1, CV_32F);
+  Mat w_full_window(averaging_radius*2 + 1, averaging_radius*2 + 1, CV_32F);
+  Mat wd_full_window(averaging_radius*2 + 1, averaging_radius*2 + 1, CV_32F);
+  float w_full_window_sum;
+
+  Mat prev_extended;
+  copyMakeBorder(prev, prev_extended, 
+                 averaging_radius, averaging_radius, averaging_radius, averaging_radius,
+                 BORDER_DEFAULT);
+
+  wd(wd_full_window, averaging_radius, averaging_radius, averaging_radius, averaging_radius, sigma_dist);
 
   for (int r0 = 0; r0 < rows; ++r0) {
     for (int c0 = 0; c0 < cols; ++c0) {
-      int u0 = floor(flow.u.at<double>(r0, c0) + 0.5);
-      int v0 = floor(flow.v.at<double>(r0, c0) + 0.5);
+      Vec2f flow_at_point = flow.at<Vec2f>(r0, c0);
+      int u0 = floor(flow_at_point[0] + 0.5);
+      if (r0 + u0 < 0) { u0 = -r0; }
+      if (r0 + u0 >= rows) { u0 = rows - 1 - r0; }
+      int v0 = floor(flow_at_point[1] + 0.5);
+      if (c0 + v0 < 0) { v0 = -c0; }
+      if (c0 + v0 >= cols) { v0 = cols - 1 - c0; }
 
       const int min_row_shift = -min(r0 + u0, max_flow);
       const int max_row_shift = min(rows - 1 - (r0 + u0), max_flow);
       const int min_col_shift = -min(c0 + v0, max_flow);
       const int max_col_shift = min(cols - 1 - (c0 + v0), max_flow);
 
-      double min_cost = DBL_MAX, best_u = u0, best_v = v0;
-
-      Mat w_full_window;
-      double w_full_window_sum;
-      Mat diff_storage;
-
-      if (r0 - averaging_radius >= 0 && 
-          r0 + averaging_radius < rows &&
-          c0 - averaging_radius >= 0 &&
-          c0 + averaging_radius < cols &&
-          mask.at<uchar>(r0, c0)) {
-          w_full_window = wd(averaging_radius, 
-                             averaging_radius, 
-                             averaging_radius, 
-                             averaging_radius, 
-                             sigma_dist).mul(
-                          wc(prev, r0, c0, 
-                             averaging_radius, 
-                             averaging_radius, 
-                             averaging_radius, 
-                             averaging_radius, 
-                             sigma_color));
+      float min_cost = DBL_MAX, best_u = u0, best_v = v0;
 
+      if (mask.at<uchar>(r0, c0)) {
+          wc(prev_extended, w_full_window, r0 + averaging_radius, c0 + averaging_radius,
+             averaging_radius, averaging_radius, averaging_radius, averaging_radius, sigma_color);
+          multiply(w_full_window, wd_full_window, w_full_window);
           w_full_window_sum = sum(w_full_window)[0];
-          diff_storage = Mat::zeros(averaging_radius*2 + 1, averaging_radius*2 + 1, CV_64F);
       }
 
       bool first_flow_iteration = true;
-      double sum_e, min_e;
+      float sum_e, min_e;
 
       for (int u = min_row_shift; u <= max_row_shift; ++u) {
         for (int v = min_col_shift; v <= max_col_shift; ++v) {
-          double e = dist(prev.at<Vec3b>(r0, c0), next.at<Vec3b>(r0 + u0 + u, c0 + v0 + v));
+          float e = dist(prev.at<Vec3b>(r0, c0), next.at<Vec3b>(r0 + u0 + u, c0 + v0 + v));
           if (first_flow_iteration) {
             sum_e = e;
             min_e = e;
@@ -269,10 +230,11 @@ static void calcOpticalFlowSingleScaleSF(const Mat& prev,
                                      r0 + u0 + u + window_bottom_shift + 1);
           const Range next_col_range(c0 + v0 + v - window_left_shift, 
                                      c0 + v0 + v + window_right_shift + 1); 
-          
+         
+                     
           Mat diff2;
           Mat w;
-          double w_sum;
+          float w_sum;
           if (window_top_shift == averaging_radius &&
               window_bottom_shift == averaging_radius &&
               window_left_shift == averaging_radius &&
@@ -280,21 +242,23 @@ static void calcOpticalFlowSingleScaleSF(const Mat& prev,
             w = w_full_window;
             w_sum = w_full_window_sum;
             diff2 = diff_storage;
-
             dist(prev(prev_row_range, prev_col_range), next(next_row_range, next_col_range), diff2);
           } else {
-            diff2 = Mat::zeros(window_bottom_shift + window_top_shift + 1,  
-                                   window_right_shift + window_left_shift + 1, CV_64F);
+            diff2 = diff_storage(Range(averaging_radius - window_top_shift, 
+                                       averaging_radius + 1 + window_bottom_shift),
+                                 Range(averaging_radius - window_left_shift,
+                                       averaging_radius + 1 + window_right_shift));
             
             dist(prev(prev_row_range, prev_col_range), next(next_row_range, next_col_range), diff2);
-
-            w = wd(window_top_shift, window_bottom_shift, window_left_shift, window_right_shift, sigma_dist).mul( 
-                wc(prev, r0, c0, window_top_shift, window_bottom_shift, window_left_shift, window_right_shift, sigma_color));
+            w = w_full_window(Range(averaging_radius - window_top_shift, 
+                                    averaging_radius + 1 + window_bottom_shift),
+                              Range(averaging_radius - window_left_shift,
+                                    averaging_radius + 1 + window_right_shift));
             w_sum = sum(w)[0];
           }
           multiply(diff2, w, diff2);
       
-          const double cost = sum(diff2)[0] / w_sum;
+          const float cost = sum(diff2)[0] / w_sum;
           if (cost < min_cost) {
             min_cost = cost;
             best_u = u + u0;
@@ -302,72 +266,37 @@ static void calcOpticalFlowSingleScaleSF(const Mat& prev,
           }
         }
       }
-      int square = (max_row_shift - min_row_shift + 1) *
-                   (max_col_shift - min_col_shift + 1);
-      confidence.at<double>(r0, c0) = (square == 0) ? 0
-                                                   : sum_e / square - min_e;
+      int windows_square = (max_row_shift - min_row_shift + 1) *
+                           (max_col_shift - min_col_shift + 1);
+      confidence.at<float>(r0, c0) = (windows_square == 0) ? 0
+                                                           : sum_e / windows_square - min_e;
+      CV_Assert(confidence.at<float>(r0, c0) >= 0); // TODO: remove it after testing
       if (mask.at<uchar>(r0, c0)) {
-        flow.u.at<double>(r0, c0) = best_u;
-        flow.v.at<double>(r0, c0) = best_v;
+        flow.at<Vec2f>(r0, c0) = Vec2f(best_u, best_v);
       }
     }
   }
 }
 
-static Flow upscaleOpticalFlow(int new_rows, 
+static Mat upscaleOpticalFlow(int new_rows, 
                                int new_cols,
                                const Mat& image,
                                const Mat& confidence,
-                               const Flow& flow,
+                               Mat& flow,
                                int averaging_radius,
-                               double sigma_dist,
-                               double sigma_color) {
-  const int rows = image.rows;
-  const int cols = image.cols;
-  Flow new_flow(new_rows, new_cols);
-  for (int r = 0; r < rows; ++r) {
-    for (int c = 0; c < cols; ++c) {
-      const int window_top_shift = min(r, averaging_radius);
-      const int window_bottom_shift = min(rows - 1 - r, averaging_radius);
-      const int window_left_shift = min(c, averaging_radius);
-      const int window_right_shift = min(cols - 1 - c, averaging_radius);
-      
-      const Range row_range(r - window_top_shift, r + window_bottom_shift + 1);
-      const Range col_range(c - window_left_shift, c + window_right_shift + 1);
-      
-      const Mat w = confidence(row_range, col_range).mul(
-        wd(window_top_shift, window_bottom_shift, window_left_shift, window_right_shift, sigma_dist)).mul( 
-        wc(image, r, c, window_top_shift, window_bottom_shift, window_left_shift, window_right_shift, sigma_color));
-
-      const double w_sum = sum(w)[0];
-      double new_u, new_v;
-      if (fabs(w_sum) < 1e-9) {
-        new_u = flow.u.at<double>(r, c);
-        new_v = flow.v.at<double>(r, c);
-      } else {
-        new_u = sum(flow.u(row_range, col_range).mul(w))[0] / w_sum;
-        new_v = sum(flow.v(row_range, col_range).mul(w))[0] / w_sum;
-      }
-      
-      for (int dr = 0; dr <= 1; ++dr) {
-        int nr = 2*r + dr;
-        for (int dc = 0; dc <= 1; ++dc) {
-          int nc = 2*c + dc;
-          if (nr < new_rows && nc < new_cols) {
-            new_flow.u.at<double>(nr, nc) = 2 * new_u;
-            new_flow.v.at<double>(nr, nc) = 2 * new_v;
-          }
-        }
-      }
-    }
-  }
+                               float sigma_dist,
+                               float sigma_color) {
+  crossBilateralFilter(flow, image, confidence, flow, averaging_radius, sigma_color, sigma_dist, false);
+  Mat new_flow;
+  resize(flow, new_flow, Size(new_cols, new_rows), 0, 0, INTER_NEAREST);
+  new_flow *= 2;
   return new_flow;
 }
 
-static Mat calcIrregularityMat(const Flow& flow, int radius) {
-  const int rows = flow.u.rows;
-  const int cols = flow.v.cols;
-  Mat irregularity = Mat::zeros(rows, cols, CV_64F);
+static Mat calcIrregularityMat(const Mat& flow, int radius) {
+  const int rows = flow.rows;
+  const int cols = flow.cols;
+  Mat irregularity(rows, cols, CV_32F);
   for (int r = 0; r < rows; ++r) {
     const int start_row = max(0, r - radius);
     const int end_row = min(rows - 1, r + radius);
@@ -376,10 +305,9 @@ static Mat calcIrregularityMat(const Flow& flow, int radius) {
       const int end_col = min(cols - 1, c + radius);
       for (int dr = start_row; dr <= end_row; ++dr) {
         for (int dc = start_col; dc <= end_col; ++dc) {
-          const double diff = dist(flow.u.at<double>(r, c), flow.v.at<double>(r, c), 
-                                   flow.u.at<double>(dr, dc), flow.v.at<double>(dr, dc));
-          if (diff > irregularity.at<double>(r, c)) {
-            irregularity.at<double>(r, c) = diff;
+          const float diff = dist(flow.at<Vec2f>(r, c), flow.at<Vec2f>(dr, dc));
+          if (diff > irregularity.at<float>(r, c)) {
+            irregularity.at<float>(r, c) = diff;
           }
         }
       }
@@ -388,7 +316,7 @@ static Mat calcIrregularityMat(const Flow& flow, int radius) {
   return irregularity;
 }
 
-static void selectPointsToRecalcFlow(const Flow& flow, 
+static void selectPointsToRecalcFlow(const Mat& flow, 
                                      int irregularity_metric_radius,
                                      int speed_up_thr,
                                      int curr_rows,
@@ -396,11 +324,10 @@ static void selectPointsToRecalcFlow(const Flow& flow,
                                      const Mat& prev_speed_up,
                                      Mat& speed_up,
                                      Mat& mask) {
-  const int prev_rows = flow.u.rows;
-  const int prev_cols = flow.v.cols;
+  const int prev_rows = flow.rows;
+  const int prev_cols = flow.cols;
 
-  Mat is_flow_regular = calcIrregularityMat(flow, 
-                                                irregularity_metric_radius)
+  Mat is_flow_regular = calcIrregularityMat(flow, irregularity_metric_radius)
                               < speed_up_thr;
   Mat done = Mat::zeros(prev_rows, prev_cols, CV_8U);
   speed_up = Mat::zeros(curr_rows, curr_cols, CV_8U);
@@ -470,28 +397,28 @@ static void selectPointsToRecalcFlow(const Flow& flow,
   }
 }
 
-static inline double extrapolateValueInRect(int height, int width,
-                                            double v11, double v12,
-                                            double v21, double v22,
+static inline float extrapolateValueInRect(int height, int width,
+                                            float v11, float v12,
+                                            float v21, float v22,
                                             int r, int c) {
   if (r == 0 && c == 0) { return v11;}
   if (r == 0 && c == width) { return v12;}
   if (r == height && c == 0) { return v21;}
   if (r == height && c == width) { return v22;}
   
-  double qr = double(r) / height;
-  double pr = 1.0 - qr;
-  double qc = double(c) / width;
-  double pc = 1.0 - qc;
+  float qr = float(r) / height;
+  float pr = 1.0 - qr;
+  float qc = float(c) / width;
+  float pc = 1.0 - qc;
 
   return v11*pr*pc + v12*pr*qc + v21*qr*pc + v22*qc*qr; 
 }
                                               
-static void extrapolateFlow(Flow& flow,
+static void extrapolateFlow(Mat& flow,
                             const Mat& speed_up) {
-  const int rows = flow.u.rows;
-  const int cols = flow.u.cols;
-  Mat done = Mat::zeros(rows, cols, CV_8U);
+  const int rows = flow.rows;
+  const int cols = flow.cols;
+  Mat done(rows, cols, CV_8U);
   for (int r = 0; r < rows; ++r) {
     for (int c = 0; c < cols; ++c) {
       if (!done.at<uchar>(r, c) && speed_up.at<uchar>(r, c) > 1) {
@@ -506,21 +433,22 @@ static void extrapolateFlow(Flow& flow,
         for (int rr = top; rr <= bottom; ++rr) {
           for (int cc = left; cc <= right; ++cc) {
             done.at<uchar>(rr, cc) = 1;
-            flow.u.at<double>(rr, cc) = extrapolateValueInRect(
-                                          height, width, 
-                                          flow.u.at<double>(top, left),
-                                          flow.u.at<double>(top, right),
-                                          flow.u.at<double>(bottom, left),
-                                          flow.u.at<double>(bottom, right),
-                                          rr-top, cc-left); 
-
-            flow.v.at<double>(rr, cc) = extrapolateValueInRect(
-                                          height, width, 
-                                          flow.v.at<double>(top, left),
-                                          flow.v.at<double>(top, right),
-                                          flow.v.at<double>(bottom, left),
-                                          flow.v.at<double>(bottom, right),
-                                          rr-top, cc-left); 
+            Vec2f flow_at_point;
+            Vec2f top_left = flow.at<Vec2f>(top, left);
+            Vec2f top_right = flow.at<Vec2f>(top, right);
+            Vec2f bottom_left = flow.at<Vec2f>(bottom, left);
+            Vec2f bottom_right = flow.at<Vec2f>(bottom, right);
+
+            flow_at_point[0] = extrapolateValueInRect(height, width, 
+                                                      top_left[0], top_right[0],
+                                                      bottom_left[0], bottom_right[0],
+                                                      rr-top, cc-left); 
+
+            flow_at_point[1] = extrapolateValueInRect(height, width, 
+                                                      top_left[1], top_right[1],
+                                                      bottom_left[1], bottom_right[1],
+                                                      rr-top, cc-left); 
+            flow.at<Vec2f>(rr, cc) = flow_at_point;
           }
         }
       }
@@ -545,8 +473,9 @@ static void buildPyramidWithResizeMethod(Mat& src,
   }
 }
 
-static Flow calcOpticalFlowSF(Mat& from, 
+void calcOpticalFlowSF(Mat& from, 
                        Mat& to, 
+                       Mat& resulted_flow,
                        int layers,
                        int averaging_block_size, 
                        int max_flow,
@@ -565,8 +494,6 @@ static Flow calcOpticalFlowSF(Mat& from,
 
   buildPyramidWithResizeMethod(from, pyr_from_images, layers - 1, INTER_CUBIC);
   buildPyramidWithResizeMethod(to, pyr_to_images, layers - 1, INTER_CUBIC);
-//  buildPyramid(from, pyr_from_images, layers - 1, BORDER_WRAP);
-//  buildPyramid(to, pyr_to_images, layers - 1, BORDER_WRAP);
 
   if ((int)pyr_from_images.size() != layers) {
       exit(1);
@@ -582,8 +509,8 @@ static Flow calcOpticalFlowSF(Mat& from,
   Mat mask = Mat::ones(first_from_image.rows, first_from_image.cols, CV_8U);
   Mat mask_inv = Mat::ones(first_from_image.rows, first_from_image.cols, CV_8U);
 
-  Flow flow(first_from_image.rows, first_from_image.cols);
-  Flow flow_inv(first_to_image.rows, first_to_image.cols);
+  Mat flow(first_from_image.rows, first_from_image.cols, CV_32FC2);
+  Mat flow_inv(first_to_image.rows, first_to_image.cols, CV_32FC2);
 
   Mat confidence;
   Mat confidence_inv;
@@ -641,8 +568,6 @@ static Flow calcOpticalFlowSF(Mat& from,
                              new_speed_up,
                              mask);
 
-    int points_to_recalculate = sum(mask)[0] / MASK_TRUE_VALUE;
-
     selectPointsToRecalcFlow(flow_inv,
                              averaging_block_size,
                              speed_up_thr,
@@ -652,8 +577,6 @@ static Flow calcOpticalFlowSF(Mat& from,
                              new_speed_up_inv,
                              mask_inv);
 
-    points_to_recalculate = sum(mask_inv)[0] / MASK_TRUE_VALUE;
-
     speed_up = new_speed_up;
     speed_up_inv = new_speed_up_inv;
 
@@ -702,55 +625,14 @@ static Flow calcOpticalFlowSF(Mat& from,
     removeOcclusions(flow_inv, flow, occ_thr, confidence_inv);
   }
 
-  WeightedCrossBilateralFilter filter_postprocess(pyr_from_images[0], 
-                                                  postprocess_window,
-                                                  sigma_dist_fix,
-                                                  sigma_color_fix);
+  crossBilateralFilter(flow, pyr_from_images[0], confidence, flow, 
+                       postprocess_window, sigma_color_fix, sigma_dist_fix);
 
-  flow.u = filter_postprocess.apply(flow.u, confidence);
-  flow.v = filter_postprocess.apply(flow.v, confidence);
-
-  Mat blured_u, blured_v;
-  GaussianBlur(flow.u, blured_u, Size(3, 3), 5);
-  GaussianBlur(flow.v, blured_v, Size(3, 3), 5);
-
-  return Flow(blured_v, blured_u);
-}
-
-void calcOpticalFlowSF(Mat& from, 
-                       Mat& to, 
-                       Mat& flowX, 
-                       Mat& flowY,
-                       int layers,
-                       int averaging_block_size, 
-                       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) {
+  GaussianBlur(flow, flow, Size(3, 3), 5);
 
-  Flow flow = calcOpticalFlowSF(from, to, 
-                                layers,
-                                averaging_block_size,
-                                max_flow,
-                                sigma_dist,
-                                sigma_color,
-                                postprocess_window,
-                                sigma_dist_fix,
-                                sigma_color_fix,
-                                occ_thr,
-                                upscale_averaging_radius,
-                                upscale_sigma_dist,
-                                upscale_sigma_color,
-                                speed_up_thr);
-  flowX = flow.u;
-  flowY = flow.v;
+  resulted_flow = Mat(flow.size(), CV_32FC2);
+  int from_to[] = {0,1 , 1,0};
+  mixChannels(&flow, 1, &resulted_flow, 1, from_to, 2);
 }
 
 }
index 55052fd..c4aa023 100644 (file)
@@ -52,32 +52,23 @@ using namespace std;
 
 namespace cv {
 
-struct Flow {
-  Mat u, v;
-
-  Flow() {;}
-
-  Flow(Mat& _u, Mat& _v)
-    : u(_u), v(_v) {;}
-  
-  Flow(int rows, int cols) {
-    u = Mat::zeros(rows, cols, CV_64F);
-    v = Mat::zeros(rows, cols, CV_64F);
-  }
-};
-
-inline static double dist(const Vec3b& p1, const Vec3b& p2) {
+inline static float dist(const Vec3b& p1, const Vec3b& p2) {
   return (p1[0] - p2[0]) * (p1[0] - p2[0]) +
          (p1[1] - p2[1]) * (p1[1] - p2[1]) +
          (p1[2] - p2[2]) * (p1[2] - p2[2]);
 }
 
-inline static double dist(const Point2f& p1, const Point2f& p2) {
+inline static float dist(const Vec2f& p1, const Vec2f& p2) {
+  return (p1[0] - p2[0]) * (p1[0] - p2[0]) +
+         (p1[1] - p2[1]) * (p1[1] - p2[1]);
+}
+
+inline static float dist(const Point2f& p1, const Point2f& p2) {
   return (p1.x - p2.x) * (p1.x - p2.x) +
          (p1.y - p2.y) * (p1.y - p2.y);
 }
 
-inline static double dist(double x1, double y1, double x2, double y2) {
+inline static float dist(float x1, float y1, float x2, float y2) {
   return (x1 - x2) * (x1 - x2) +
          (y1 - y2) * (y1 - y2);
 }
@@ -92,34 +83,6 @@ inline static T min(T t1, T t2, T t3) {
   return (t1 <= t2 && t1 <= t3) ? t1 : min(t2, t3);
 }
 
-template<class T>
-vector<vector<T> > build(int n, int m) {
-  vector<vector<T> > res(n);
-  for (int i = 0; i < n; ++i) {
-    res[i].resize(m, 0);
-  }
-  return res;
-}
-
-class WeightedCrossBilateralFilter {
-public:
-  WeightedCrossBilateralFilter(const Mat& _image,
-                       int _windowSize,
-                       double _sigmaDist,
-                       double _sigmaColor);
-
-  Mat apply(Mat& matrix, Mat& weights);
-
-private:
-  double convolution(Mat& matrix, int row, int col, Mat& weights);
-
-  Mat image;
-  int windowSize;
-  double sigmaDist, sigmaColor;
-
-  vector<double> expDist;
-  vector<vector<vector<vector<double> > > > wc;
-};
 }
 
 #endif
index 186ba8f..050d595 100644 (file)
@@ -58,66 +58,67 @@ protected:
 
 CV_SimpleFlowTest::CV_SimpleFlowTest() {}
 
-static void readOpticalFlowFromFile(FILE* file, cv::Mat& flowX, cv::Mat& flowY) {
+static bool readOpticalFlowFromFile(FILE* file, cv::Mat& flow) {
   char header[5];
   if (fread(header, 1, 4, file) < 4 && (string)header != "PIEH") {
-    return;
+    return false;
   }
 
   int cols, rows;
   if (fread(&cols, sizeof(int), 1, file) != 1||
       fread(&rows, sizeof(int), 1, file) != 1) {
-    return;
+    return false;
   }
 
-  flowX = cv::Mat::zeros(rows, cols, CV_64F);
-  flowY = cv::Mat::zeros(rows, cols, CV_64F);
+  flow = cv::Mat::zeros(rows, cols, CV_32FC2);
 
   for (int i = 0; i < rows; ++i) {
     for (int j = 0; j < cols; ++j) {
-      float uPoint, vPoint;
-      if (fread(&uPoint, sizeof(float), 1, file) != 1 ||
-          fread(&vPoint, sizeof(float), 1, file) != 1) {
-        flowX.release();
-        flowY.release();
-        return;
+      cv::Vec2f flow_at_point;
+      if (fread(&(flow_at_point[0]), sizeof(float), 1, file) != 1 ||
+          fread(&(flow_at_point[1]), sizeof(float), 1, file) != 1) {
+        return false;
       }
-  
-      flowX.at<double>(i, j) = uPoint;
-      flowY.at<double>(i, j) = vPoint;
+      flow.at<cv::Vec2f>(i, j) = flow_at_point;
     }
   }
+
+  return true;
 }
 
-static bool isFlowCorrect(double u) {
+static bool isFlowCorrect(float u) {
   return !isnan(u) && (fabs(u) < 1e9);
 }
 
-static double calc_rmse(cv::Mat flow1X, cv::Mat flow1Y, cv::Mat flow2X, cv::Mat flow2Y) {
-  long double sum;
+static float calc_rmse(cv::Mat flow1, cv::Mat flow2) {
+  float sum;
   int counter = 0;
-  const int rows = flow1X.rows;
-  const int cols = flow1X.cols;
+  const int rows = flow1.rows;
+  const int cols = flow1.cols;
 
   for (int y = 0; y < rows; ++y) {
     for (int x = 0; x < cols; ++x) {
-      double u1 = flow1X.at<double>(y, x);
-      double v1 = flow1Y.at<double>(y, x);
-      double u2 = flow2X.at<double>(y, x);
-      double v2 = flow2Y.at<double>(y, x);
+      cv::Vec2f flow1_at_point = flow1.at<cv::Vec2f>(y, x);
+      cv::Vec2f flow2_at_point = flow2.at<cv::Vec2f>(y, x);
+
+      float u1 = flow1_at_point[0];
+      float v1 = flow1_at_point[1];
+      float u2 = flow2_at_point[0];
+      float v2 = flow2_at_point[1];
+
       if (isFlowCorrect(u1) && isFlowCorrect(u2) && isFlowCorrect(v1) && isFlowCorrect(v2)) {
         sum += (u1-u2)*(u1-u2) + (v1-v2)*(v1-v2);
         counter++;
       }
     }
   }
-  return sqrt((double)sum / (1e-9 + counter));
+  return sqrt(sum / (1e-9 + counter));
 }
 
 void CV_SimpleFlowTest::run(int) {
     int code = cvtest::TS::OK;
     
-    const double MAX_RMSE = 0.6;
+    const float MAX_RMSE = 0.6;
     const string frame1_path = ts->get_data_path() + "optflow/RubberWhale1.png";
     const string frame2_path = ts->get_data_path() + "optflow/RubberWhale2.png";
     const string gt_flow_path = ts->get_data_path() + "optflow/RubberWhale.flo";
@@ -151,7 +152,7 @@ void CV_SimpleFlowTest::run(int) {
       return;
     }
 
-    cv::Mat flowX_gt, flowY_gt;
+    cv::Mat flow_gt;
 
     FILE* gt_flow_file = fopen(gt_flow_path.c_str(), "rb");
     if (gt_flow_file == NULL) {
@@ -160,8 +161,8 @@ void CV_SimpleFlowTest::run(int) {
       ts->set_failed_test_info(cvtest::TS::FAIL_MISSING_TEST_DATA);
       return;
     }
-    readOpticalFlowFromFile(gt_flow_file, flowX_gt, flowY_gt);
-    if (flowX_gt.empty() || flowY_gt.empty()) {
+
+    if (!readOpticalFlowFromFile(gt_flow_file, flow_gt)) {
       ts->printf(cvtest::TS::LOG, "error while reading flow data from file %s",
                  gt_flow_path.c_str());
       ts->set_failed_test_info(cvtest::TS::FAIL_MISSING_TEST_DATA);
@@ -169,12 +170,12 @@ void CV_SimpleFlowTest::run(int) {
     }
     fclose(gt_flow_file);
 
-    cv::Mat flowX, flowY;
+    cv::Mat flow;
     cv::calcOpticalFlowSF(frame1, frame2, 
-                          flowX, flowY,
+                          flow,
                           3, 4, 2, 4.1, 25.5, 18, 55.0, 25.5, 0.35, 18, 55.0, 25.5, 10);
 
-    double rmse = calc_rmse(flowX_gt, flowY_gt, flowX, flowY);
+    float rmse = calc_rmse(flow_gt, flow);
     
     ts->printf(cvtest::TS::LOG, "Optical flow estimation RMSE for SimpleFlow algorithm : %lf\n",
                rmse);
index 6a195fe..332df78 100644 (file)
 using namespace cv;
 using namespace std;
 
+#define APP_NAME "simpleflow_demo : "
+
 static void help()
 {
-    // print a welcome message, and the OpenCV version
-    printf("This is a demo of SimpleFlow optical flow algorithm,\n"
-           "Using OpenCV version %s\n\n", CV_VERSION);
-
-    printf("Usage: simpleflow_demo frame1 frame2 output_flow"
-           "\nApplication will write estimated flow "
-           "\nbetween 'frame1' and 'frame2' in binary format"
-           "\ninto file 'output_flow'"
-           "\nThen one can use code from http://vision.middlebury.edu/flow/data/"
-           "\nto convert flow in binary file to image\n");
+  // print a welcome message, and the OpenCV version
+  printf("This is a demo of SimpleFlow optical flow algorithm,\n"
+         "Using OpenCV version %s\n\n", CV_VERSION);
+
+  printf("Usage: simpleflow_demo frame1 frame2 output_flow"
+         "\nApplication will write estimated flow "
+         "\nbetween 'frame1' and 'frame2' in binary format"
+         "\ninto file 'output_flow'"
+         "\nThen one can use code from http://vision.middlebury.edu/flow/data/"
+         "\nto convert flow in binary file to image\n");
 }
 
 // binary file format for flow data specified here:
 // http://vision.middlebury.edu/flow/data/
-static void writeOpticalFlowToFile(const Mat& u, const Mat& v, FILE* file) {
-  int cols = u.cols;
-  int rows = u.rows;
+static void writeOpticalFlowToFile(const Mat& flow, FILE* file) {
+  int cols = flow.cols;
+  int rows = flow.rows;
   
   fprintf(file, "PIEH");
-   
   if (fwrite(&cols, sizeof(int), 1, file) != 1 ||
       fwrite(&rows, sizeof(int), 1, file) != 1) {
-    fprintf(stderr, "writeOpticalFlowToFile : problem writing header\n");
+    printf(APP_NAME "writeOpticalFlowToFile : problem writing header\n");
     exit(1);
   }
 
-  for (int i= 0; i < u.rows; ++i) {
-    for (int j = 0; j < u.cols; ++j) {
-      float uPoint = u.at<double>(i, j);
-      float vPoint = v.at<double>(i, j);
+  for (int i= 0; i < rows; ++i) {
+    for (int j = 0; j < cols; ++j) {
+      Vec2f flow_at_point = flow.at<Vec2f>(i, j);
 
-      if (fwrite(&uPoint, sizeof(float), 1, file) != 1 ||
-          fwrite(&vPoint, sizeof(float), 1, file) != 1) {
-        fprintf(stderr, "writeOpticalFlowToFile : problem writing data\n");
+      if (fwrite(&(flow_at_point[0]), sizeof(float), 1, file) != 1 ||
+          fwrite(&(flow_at_point[1]), sizeof(float), 1, file) != 1) {
+        printf(APP_NAME "writeOpticalFlowToFile : problem writing data\n");
         exit(1);
       }
     }
   }
 }
-int main(int argc, char** argv) {
-    help();
 
-    if (argc < 4) {
-      fprintf(stderr, "Wrong number of command line arguments : %d (expected %d)\n", argc, 4);
-      exit(1);
-    }
-    
-    Mat frame1 = imread(argv[1]);
-    Mat frame2 = imread(argv[2]);
+static void run(int argc, char** argv) {
+  if (argc < 3) {
+    printf(APP_NAME "Wrong number of command line arguments for mode `run`: %d (expected %d)\n", 
+           argc, 3);
+    exit(1);
+  }
 
-    if (frame1.empty() || frame2.empty()) {
-      fprintf(stderr, "simpleflow_demo : Images cannot be read\n");
-      exit(1);
-    }
+  Mat frame1 = imread(argv[0]);
+  Mat frame2 = imread(argv[1]);
 
-    if (frame1.rows != frame2.rows && frame1.cols != frame2.cols) {
-      fprintf(stderr, "simpleflow_demo : Images should be of equal sizes\n");
-      exit(1);
-    }
+  if (frame1.empty()) {
+    printf(APP_NAME "Image #1 : %s cannot be read\n", argv[0]);
+    exit(1);
+  }
 
-    if (frame1.type() != 16 || frame2.type() != 16) {
-      fprintf(stderr, "simpleflow_demo : Images should be of equal type CV_8UC3\n");
-      exit(1);
-    }
+  if (frame2.empty()) {
+    printf(APP_NAME "Image #2 : %s cannot be read\n", argv[1]);
+    exit(1);
+  }
+
+  if (frame1.rows != frame2.rows && frame1.cols != frame2.cols) {
+    printf(APP_NAME "Images should be of equal sizes\n");
+    exit(1);
+  }
 
-    printf("simpleflow_demo : Read two images of size [rows = %d, cols = %d]\n", 
-           frame1.rows, frame1.cols);
+  if (frame1.type() != 16 || frame2.type() != 16) {
+    printf(APP_NAME "Images should be of equal type CV_8UC3\n");
+    exit(1);
+  }
+
+  printf(APP_NAME "Read two images of size [rows = %d, cols = %d]\n", 
+         frame1.rows, frame1.cols);
 
-    Mat flowX, flowY;
+  Mat flow;
 
-    calcOpticalFlowSF(frame1, frame2, 
-                      flowX, flowY,
-                      3, 2, 4, 4.1, 25.5, 18, 55.0, 25.5, 0.35, 18, 55.0, 25.5, 10);
+  float start = getTickCount();
+  calcOpticalFlowSF(frame1, frame2, 
+                    flow,
+                    3, 2, 4, 4.1, 25.5, 18, 55.0, 25.5, 0.35, 18, 55.0, 25.5, 10);
+  printf(APP_NAME "calcOpticalFlowSF : %lf sec\n", (getTickCount() - start) / getTickFrequency());
 
-  FILE* file = fopen(argv[3], "wb");
+  FILE* file = fopen(argv[2], "wb");
   if (file == NULL) {
-    fprintf(stderr, "simpleflow_demo : Unable to open file '%s' for writing\n", argv[3]); 
+    printf(APP_NAME "Unable to open file '%s' for writing\n", argv[2]); 
     exit(1);
   }
-  printf("simpleflow_demo : Writing to file\n");
-  writeOpticalFlowToFile(flowX, flowY, file);
+  printf(APP_NAME "Writing to file\n");
+  writeOpticalFlowToFile(flow, file);
   fclose(file);
+}
+
+static bool readOpticalFlowFromFile(FILE* file, Mat& flow) {
+  char header[5];
+  if (fread(header, 1, 4, file) < 4 && (string)header != "PIEH") {
+    return false;
+  }
+
+  int cols, rows;
+  if (fread(&cols, sizeof(int), 1, file) != 1||
+      fread(&rows, sizeof(int), 1, file) != 1) {
+    return false;
+  }
+
+  flow = Mat::zeros(rows, cols, CV_32FC2);
+
+  for (int i = 0; i < rows; ++i) {
+    for (int j = 0; j < cols; ++j) {
+      Vec2f flow_at_point;
+      if (fread(&(flow_at_point[0]), sizeof(float), 1, file) != 1 ||
+          fread(&(flow_at_point[1]), sizeof(float), 1, file) != 1) {
+        return false;
+      }
+      flow.at<Vec2f>(i, j) = flow_at_point;
+    }
+  }
+
+  return true;
+}
+
+static bool isFlowCorrect(float u) {
+  return !isnan(u) && (fabs(u) < 1e9);
+}
+
+static float calc_rmse(Mat flow1, Mat flow2) {
+  float sum;
+  int counter = 0;
+  const int rows = flow1.rows;
+  const int cols = flow1.cols;
+
+  for (int y = 0; y < rows; ++y) {
+    for (int x = 0; x < cols; ++x) {
+      Vec2f flow1_at_point = flow1.at<Vec2f>(y, x);
+      Vec2f flow2_at_point = flow2.at<Vec2f>(y, x);
+
+      float u1 = flow1_at_point[0];
+      float v1 = flow1_at_point[1];
+      float u2 = flow2_at_point[0];
+      float v2 = flow2_at_point[1];
+
+      if (isFlowCorrect(u1) && isFlowCorrect(u2) && isFlowCorrect(v1) && isFlowCorrect(v2)) {
+        sum += (u1-u2)*(u1-u2) + (v1-v2)*(v1-v2);
+        counter++;
+      }
+    }
+  }
+  return sqrt(sum / (1e-9 + counter));
+}
+
+static void eval(int argc, char** argv) {
+  if (argc < 2) {
+    printf(APP_NAME "Wrong number of command line arguments for mode `eval` : %d (expected %d)\n", 
+           argc, 2);
+    exit(1);
+  }
+  
+  Mat flow1, flow2;
+
+  FILE* flow_file_1 = fopen(argv[0], "rb");
+  if (flow_file_1 == NULL) {
+    printf(APP_NAME "Cannot open file with first flow : %s\n", argv[0]);
+    exit(1);
+  }
+  if (!readOpticalFlowFromFile(flow_file_1, flow1)) {
+    printf(APP_NAME "Cannot read flow data from file %s\n", argv[0]);
+    exit(1);
+  }
+  fclose(flow_file_1);
+
+  FILE* flow_file_2 = fopen(argv[1], "rb");
+  if (flow_file_2 == NULL) {
+    printf(APP_NAME "Cannot open file with first flow : %s\n", argv[1]);
+    exit(1);
+  }
+  if (!readOpticalFlowFromFile(flow_file_2, flow2)) {
+    printf(APP_NAME "Cannot read flow data from file %s\n", argv[1]);
+    exit(1);
+  }
+  fclose(flow_file_2);
+  
+  float rmse = calc_rmse(flow1, flow2);
+  printf("%lf\n", rmse);
+}
+
+int main(int argc, char** argv) {
+  if (argc < 2) {
+    printf(APP_NAME "Mode is not specified\n");
+    help();
+    exit(1);
+  }
+  string mode = (string)argv[1];
+  int new_argc = argc - 2;
+  char** new_argv = &argv[2];
+
+  if ("run" == mode) {
+    run(new_argc, new_argv);
+  } else if ("eval" == mode) {
+    eval(new_argc, new_argv);
+  } else if ("help" == mode) 
+    help();
+  else {
+    printf(APP_NAME "Unknown mode : %s\n", argv[1]);   
+    help();
+  }
+  
   return 0;
 }