From 784c09d6f9f85f1aa0774fdca56387cd6817868d Mon Sep 17 00:00:00 2001 From: Yury Zemlyanskiy Date: Mon, 27 Aug 2012 14:40:57 +0400 Subject: [PATCH] Updates for SimpleFlow algorithm + 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 | 5 +- modules/video/src/simpleflow.cpp | 450 +++++++++-------------- modules/video/src/simpleflow.hpp | 53 +-- modules/video/test/test_simpleflow.cpp | 63 ++-- samples/cpp/simpleflow_demo.cpp | 229 +++++++++--- 5 files changed, 384 insertions(+), 416 deletions(-) diff --git a/modules/video/include/opencv2/video/tracking.hpp b/modules/video/include/opencv2/video/tracking.hpp index 85c1881..6800c63 100644 --- a/modules/video/include/opencv2/video/tracking.hpp +++ b/modules/video/include/opencv2/video/tracking.hpp @@ -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, diff --git a/modules/video/src/simpleflow.cpp b/modules/video/src/simpleflow.cpp index 1fda361..0cd320d 100644 --- a/modules/video/src/simpleflow.cpp +++ b/modules/video/src/simpleflow.cpp @@ -54,193 +54,154 @@ 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(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(row, col), - image.at(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(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(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(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(r, c), flow.v.at(r, c), - -flow_inv.u.at(r, c), -flow_inv.v.at(r, c)) > occ_thr) { - confidence.at(r, c) = 0; - occlusions++; + if (dist(flow.at(r, c), -flow_inv.at(r, c)) > occ_thr) { + confidence.at(r, c) = 0; + } else { + confidence.at(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(r, c) = -(dr*dr + dc*dc) * factor; + d.at(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(r0, c0); for (int dr = r0-top_shift, r = 0; dr <= r0+bottom_shift; ++dr, ++r) { + const Vec3b *row = image.ptr(dr); + float *d_row = d.ptr(r); for (int dc = c0-left_shift, c = 0; dc <= c0+right_shift; ++dc, ++c) { - d.at(r, c) = -dist(image.at(r0, c0), image.at(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(r); const Vec3b *m2_row = m2.ptr(r); - double* row = result.ptr(r); + float* row = result.ptr(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 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(row, col)[ch] = (flag && fabs(weights_sum) < 1e-9) + ? image.at(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(r0, c0) + 0.5); - int v0 = floor(flow.v.at(r0, c0) + 0.5); + Vec2f flow_at_point = flow.at(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(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(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(r0, c0), next.at(r0 + u0 + u, c0 + v0 + v)); + float e = dist(prev.at(r0, c0), next.at(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(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(r0, c0) = (windows_square == 0) ? 0 + : sum_e / windows_square - min_e; + CV_Assert(confidence.at(r0, c0) >= 0); // TODO: remove it after testing if (mask.at(r0, c0)) { - flow.u.at(r0, c0) = best_u; - flow.v.at(r0, c0) = best_v; + flow.at(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(r, c); - new_v = flow.v.at(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(nr, nc) = 2 * new_u; - new_flow.v.at(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(r, c), flow.v.at(r, c), - flow.u.at(dr, dc), flow.v.at(dr, dc)); - if (diff > irregularity.at(r, c)) { - irregularity.at(r, c) = diff; + const float diff = dist(flow.at(r, c), flow.at(dr, dc)); + if (diff > irregularity.at(r, c)) { + irregularity.at(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(r, c) && speed_up.at(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(rr, cc) = 1; - flow.u.at(rr, cc) = extrapolateValueInRect( - height, width, - flow.u.at(top, left), - flow.u.at(top, right), - flow.u.at(bottom, left), - flow.u.at(bottom, right), - rr-top, cc-left); - - flow.v.at(rr, cc) = extrapolateValueInRect( - height, width, - flow.v.at(top, left), - flow.v.at(top, right), - flow.v.at(bottom, left), - flow.v.at(bottom, right), - rr-top, cc-left); + Vec2f flow_at_point; + Vec2f top_left = flow.at(top, left); + Vec2f top_right = flow.at(top, right); + Vec2f bottom_left = flow.at(bottom, left); + Vec2f bottom_right = flow.at(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(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); } } diff --git a/modules/video/src/simpleflow.hpp b/modules/video/src/simpleflow.hpp index 55052fd..c4aa023 100644 --- a/modules/video/src/simpleflow.hpp +++ b/modules/video/src/simpleflow.hpp @@ -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 -vector > build(int n, int m) { - vector > 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 expDist; - vector > > > wc; -}; } #endif diff --git a/modules/video/test/test_simpleflow.cpp b/modules/video/test/test_simpleflow.cpp index 186ba8f..050d595 100644 --- a/modules/video/test/test_simpleflow.cpp +++ b/modules/video/test/test_simpleflow.cpp @@ -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(i, j) = uPoint; - flowY.at(i, j) = vPoint; + flow.at(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(y, x); - double v1 = flow1Y.at(y, x); - double u2 = flow2X.at(y, x); - double v2 = flow2Y.at(y, x); + cv::Vec2f flow1_at_point = flow1.at(y, x); + cv::Vec2f flow2_at_point = flow2.at(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); diff --git a/samples/cpp/simpleflow_demo.cpp b/samples/cpp/simpleflow_demo.cpp index 6a195fe..332df78 100644 --- a/samples/cpp/simpleflow_demo.cpp +++ b/samples/cpp/simpleflow_demo.cpp @@ -8,89 +8,212 @@ 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(i, j); - float vPoint = v.at(i, j); + for (int i= 0; i < rows; ++i) { + for (int j = 0; j < cols; ++j) { + Vec2f flow_at_point = flow.at(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(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(y, x); + Vec2f flow2_at_point = flow2.at(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; } -- 2.7.4