Skip to content

O
opencv_contrib
Commit 3ff3cd84 authored by Maksim Shabunin's avatar Maksim Shabunin
Browse files
Options

Merge pull request #689 from sbokov:DIS_optical_flow

parents cd5993c6 d4c3765e
Hide whitespace changes
Inline Side-by-side
Showing with 791 additions and 2 deletions
modules/optflow/doc/optflow.bib
View file @ 3ff3cd84
......@@ -37,3 +37,10 @@
year={2013},
organization={IEEE}
}
@article{Kroeger2016,
title={Fast Optical Flow using Dense Inverse Search},
author={Kroeger, Till and Timofte, Radu and Dai, Dengxin and Van Gool, Luc},
journal={arXiv preprint arXiv:1603.03590},
year={2016}
}
modules/optflow/include/opencv2/optflow.hpp
View file @ 3ff3cd84
......@@ -191,6 +191,44 @@ CV_EXPORTS_W Ptr<DenseOpticalFlow> createOptFlow_Farneback();
//! Additional interface to the SparseToDenseFlow algorithm - calcOpticalFlowSparseToDense()
CV_EXPORTS_W Ptr<DenseOpticalFlow> createOptFlow_SparseToDense();
/** @brief DIS optical flow algorithm.
This class implements the Dense Inverse Search (DIS) optical flow algorithm. More
details about the algorithm can be found at @cite Kroeger2016 .
*/
class CV_EXPORTS_W DISOpticalFlow : public DenseOpticalFlow
{
public:
/** @brief Finest level of the gaussian pyramid on which the flow is computed (zero level
corresponds to the original image resolution).The final flow is obtained by bilinear upscaling.
@see setFinestScale */
virtual int getFinestScale() const = 0;
/** @copybrief getFinestScale @see getFinestScale */
virtual void setFinestScale(int val) = 0;
/** @brief Size of an image patch for matching (in pixels)
@see setPatchSize */
virtual int getPatchSize() const = 0;
/** @copybrief getPatchSize @see getPatchSize */
virtual void setPatchSize(int val) = 0;
/** @brief Stride between neighbor patches. Must be less than patch size.
@see setPatchStride */
virtual int getPatchStride() const = 0;
/** @copybrief getPatchStride @see getPatchStride */
virtual void setPatchStride(int val) = 0;
/** @brief number of gradient descent iterations in the patch inverse search stage
@see setGradientDescentIterations */
virtual int getGradientDescentIterations() const = 0;
/** @copybrief getGradientDescentIterations @see getGradientDescentIterations */
virtual void setGradientDescentIterations(int val) = 0;
};
/** @brief Creates an instance of DISOpticalFlow
*/
CV_EXPORTS_W Ptr<DISOpticalFlow> createOptFlow_DIS();
//! @}
} //optflow
......
modules/optflow/samples/optical_flow_benchmark.py 0 → 100644
View file @ 3ff3cd84
#!/usr/bin/env python
from __future__ import print_function
import os, sys, shutil
import argparse
import json, re
from subprocess import check_output
import datetime
import matplotlib.pyplot as plt
def load_json(path):
f = open(path, "r")
data = json.load(f)
return data
def save_json(obj, path):
tmp_file = path + ".bak"
f = open(tmp_file, "w")
json.dump(obj, f, indent=2)
f.flush()
os.fsync(f.fileno())
f.close()
try:
os.rename(tmp_file, path)
except:
os.remove(path)
os.rename(tmp_file, path)
def parse_evaluation_result(input_str, i):
res = {}
res['frame_number'] = i + 1
res['error'] = {}
regex = "([A-Za-z. \\[\\].0-9]+):[ ]*([0-9]*\.[0-9]+|[0-9]+)"
for elem in re.findall(regex,input_str):
if "Time" in elem[0]:
res['time'] = float(elem[1])
elif "Average" in elem[0]:
res['error']['average'] = float(elem[1])
elif "deviation" in elem[0]:
res['error']['std'] = float(elem[1])
else:
res['error'][elem[0]] = float(elem[1])
return res
def evaluate_sequence(sequence, algorithm, dataset, executable, img_files, gt_files,
state, state_path):
if "eval_results" not in state[dataset][algorithm][-1].keys():
state[dataset][algorithm][-1]["eval_results"] = {}
elif sequence in state[dataset][algorithm][-1]["eval_results"].keys():
return
res = []
for i in range(len(img_files) - 1):
sys.stdout.write("Algorithm: %-20s Sequence: %-10s Done: [%3d/%3d]\r" %
(algorithm, sequence, i, len(img_files) - 1)),
sys.stdout.flush()
res_string = check_output([executable, img_files[i], img_files[i + 1],
algorithm, gt_files[i]])
res.append(parse_evaluation_result(res_string, i))
state[dataset][algorithm][-1]["eval_results"][sequence] = res
save_json(state, state_path)
#############################DATSET DEFINITIONS################################
def evaluate_mpi_sintel(source_dir, algorithm, evaluation_executable, state, state_path):
evaluation_result = {}
img_dir = os.path.join(source_dir, 'mpi_sintel', 'training', 'final')
gt_dir = os.path.join(source_dir, 'mpi_sintel', 'training', 'flow')
sequences = [f for f in os.listdir(img_dir)
if os.path.isdir(os.path.join(img_dir, f))]
for seq in sequences:
img_files = sorted([os.path.join(img_dir, seq, f)
for f in os.listdir(os.path.join(img_dir, seq))
if f.endswith(".png")])
gt_files = sorted([os.path.join(gt_dir, seq, f)
for f in os.listdir(os.path.join(gt_dir, seq))
if f.endswith(".flo")])
evaluation_result[seq] = evaluate_sequence(seq, algorithm, 'mpi_sintel',
evaluation_executable, img_files, gt_files, state, state_path)
return evaluation_result
def evaluate_middlebury(source_dir, algorithm, evaluation_executable, state, state_path):
evaluation_result = {}
img_dir = os.path.join(source_dir, 'middlebury', 'other-data')
gt_dir = os.path.join(source_dir, 'middlebury', 'other-gt-flow')
sequences = [f for f in os.listdir(gt_dir)
if os.path.isdir(os.path.join(gt_dir, f))]
for seq in sequences:
img_files = sorted([os.path.join(img_dir, seq, f)
for f in os.listdir(os.path.join(img_dir, seq))
if f.endswith(".png")])
gt_files = sorted([os.path.join(gt_dir, seq, f)
for f in os.listdir(os.path.join(gt_dir, seq))
if f.endswith(".flo")])
evaluation_result[seq] = evaluate_sequence(seq, algorithm, 'middlebury',
evaluation_executable, img_files, gt_files, state, state_path)
return evaluation_result
dataset_eval_functions = {
"mpi_sintel": evaluate_mpi_sintel,
"middlebury": evaluate_middlebury
}
###############################################################################
def create_dir(dir):
if not os.path.exists(dir):
os.makedirs(dir)
def parse_sequence(input_str):
if len(input_str) == 0:
return []
else:
return [o.strip() for o in input_str.split(",") if o]
def build_chart(dst_folder, state, dataset):
fig = plt.figure(figsize=(16, 10))
markers = ["o", "s", "h", "^", "D"]
marker_idx = 0
colors = ["b", "g", "r"]
color_idx = 0
for algo in state[dataset].keys():
for eval_instance in state[dataset][algo]:
name = algo + "--" + eval_instance["timestamp"]
average_time = 0.0
average_error = 0.0
num_elem = 0
for seq in eval_instance["eval_results"].keys():
for frame in eval_instance["eval_results"][seq]:
average_time += frame["time"]
average_error += frame["error"]["average"]
num_elem += 1
average_time /= num_elem
average_error /= num_elem
marker_style = colors[color_idx] + markers[marker_idx]
color_idx += 1
if color_idx >= len(colors):
color_idx = 0
marker_idx += 1
if marker_idx >= len(markers):
marker_idx = 0
plt.gca().plot([average_time], [average_error],
marker_style,
markersize=14,
label=name)
plt.gca().set_ylabel('Average Endpoint Error (EPE)', fontsize=20)
plt.gca().set_xlabel('Average Runtime (seconds per frame)', fontsize=20)
plt.gca().set_xscale("log")
plt.gca().set_title('Evaluation on ' + dataset, fontsize=20)
plt.gca().legend()
fig.savefig(os.path.join(dst_folder, "evaluation_results_" + dataset + ".png"),
bbox_inches='tight')
plt.close()
if __name__ == '__main__':
parser = argparse.ArgumentParser(
description='Optical flow benchmarking script',
formatter_class=argparse.RawDescriptionHelpFormatter)
parser.add_argument(
"bin_path",
default="./optflow-example-optical_flow_evaluation",
help="Path to the optical flow evaluation executable")
parser.add_argument(
"-a",
"--algorithms",
metavar="ALGORITHMS",
default="",
help=("Comma-separated list of optical-flow algorithms to evaluate "
"(example: -a farneback,tvl1,deepflow). Note that previously "
"evaluated algorithms are also included in the output charts"))
parser.add_argument(
"-d",
"--datasets",
metavar="DATASETS",
default="mpi_sintel",
help=("Comma-separated list of datasets for evaluation (currently only "
"'mpi_sintel' and 'middlebury' are supported)"))
parser.add_argument(
"-f",
"--dataset_folder",
metavar="DATASET_FOLDER",
default="./OF_datasets",
help=("Path to a folder containing datasets. To enable evaluation on "
"MPI Sintel dataset, please download it using the following links: "
"http://files.is.tue.mpg.de/sintel/MPI-Sintel-training_images.zip and "
"http://files.is.tue.mpg.de/sintel/MPI-Sintel-training_extras.zip and "
"unzip these archives into the 'mpi_sintel' folder. To enable evaluation "
"on the Middlebury dataset use the following links: "
"http://vision.middlebury.edu/flow/data/comp/zip/other-color-twoframes.zip, "
"http://vision.middlebury.edu/flow/data/comp/zip/other-gt-flow.zip. "
"These should be unzipped into 'middlebury' folder"))
parser.add_argument(
"-o",
"--out",
metavar="OUT_DIR",
default="./OF_evaluation_results",
help="Output directory where to store benchmark results")
parser.add_argument(
"-s",
"--state",
metavar="STATE_JSON",
default="./OF_evaluation_state.json",
help=("Path to a json file that stores the current evaluation state and "
"previous evaluation results"))
args, other_args = parser.parse_known_args()
if not os.path.isfile(args.bin_path):
print("Error: " + args.bin_path + " does not exist")
sys.exit(1)
if not os.path.exists(args.dataset_folder):
print("Error: " + args.dataset_folder + (" does not exist. Please, correctly "
"specify the -f parameter"))
sys.exit(1)
state = {}
if os.path.isfile(args.state):
state = load_json(args.state)
algorithm_list = parse_sequence(args.algorithms)
dataset_list = parse_sequence(args.datasets)
for dataset in dataset_list:
if dataset not in dataset_eval_functions.keys():
print("Error: unsupported dataset " + dataset)
sys.exit(1)
if dataset not in os.listdir(args.dataset_folder):
print("Error: " + os.path.join(args.dataset_folder, dataset) + (" does not exist. "
"Please, download the dataset and follow the naming conventions "
"(use -h for more information)"))
sys.exit(1)
for dataset in dataset_list:
if dataset not in state.keys():
state[dataset] = {}
for algorithm in algorithm_list:
if algorithm in state[dataset].keys():
last_eval_instance = state[dataset][algorithm][-1]
if "finished" not in last_eval_instance.keys():
print(("Continuing an unfinished evaluation of " +
algorithm + " started at " + last_eval_instance["timestamp"]))
else:
state[dataset][algorithm].append({"timestamp":
datetime.datetime.now().strftime("%Y-%m-%d--%H-%M")})
else:
state[dataset][algorithm] = [{"timestamp":
datetime.datetime.now().strftime("%Y-%m-%d--%H-%M")}]
save_json(state, args.state)
dataset_eval_functions[dataset](args.dataset_folder, algorithm, args.bin_path,
state, args.state)
state[dataset][algorithm][-1]["finished"] = True
save_json(state, args.state)
save_json(state, args.state)
create_dir(args.out)
for dataset in dataset_list:
build_chart(args.out, state, dataset)
modules/optflow/samples/optical_flow_evaluation.cpp
View file @ 3ff3cd84
......@@ -11,7 +11,7 @@ using namespace optflow;
const String keys = "{help h usage ? | | print this message }"
"{@image1 | | image1 }"
"{@image2 | | image2 }"
"{@algorithm | | [farneback, simpleflow, tvl1, deepflow or sparsetodenseflow] }"
"{@algorithm | | [farneback, simpleflow, tvl1, deepflow, sparsetodenseflow or DISflow] }"
"{@groundtruth | | path to the .flo file (optional), Middlebury format }"
"{m measure |endpoint| error measure - [endpoint or angular] }"
"{r region |all | region to compute stats about [all, discontinuities, untextured] }"
......@@ -229,7 +229,7 @@ int main( int argc, char** argv )
if ( i2.depth() != CV_8U )
i2.convertTo(i2, CV_8U);
if ( (method == "farneback" || method == "tvl1" || method == "deepflow") && i1.channels() == 3 )
if ( (method == "farneback" || method == "tvl1" || method == "deepflow" || method == "DISflow") && i1.channels() == 3 )
{ // 1-channel images are expected
cvtColor(i1, i1, COLOR_BGR2GRAY);
cvtColor(i2, i2, COLOR_BGR2GRAY);
......@@ -252,6 +252,8 @@ int main( int argc, char** argv )
algorithm = createOptFlow_DeepFlow();
else if ( method == "sparsetodenseflow" )
algorithm = createOptFlow_SparseToDense();
else if ( method == "DISflow" )
algorithm = createOptFlow_DIS();
else
{
printf("Wrong method!\n");
......
modules/optflow/src/dis_flow.cpp 0 → 100644
View file @ 3ff3cd84
/*M///////////////////////////////////////////////////////////////////////////////////////
//
// IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
//
// By downloading, copying, installing or using the software you agree to this license.
// If you do not agree to this license, do not download, install,
// copy or use the software.
//
//
// License Agreement
// For Open Source Computer Vision Library
//
// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
// Third party copyrights are property of their respective owners.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
// * Redistribution's of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
//
// * Redistribution's in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
//
// * The name of the copyright holders may not be used to endorse or promote products
// derived from this software without specific prior written permission.
//
// This software is provided by the copyright holders and contributors "as is" and
// any express or implied warranties, including, but not limited to, the implied
// warranties of merchantability and fitness for a particular purpose are disclaimed.
// In no event shall the Intel Corporation or contributors be liable for any direct,
// indirect, incidental, special, exemplary, or consequential damages
// (including, but not limited to, procurement of substitute goods or services;
// loss of use, data, or profits; or business interruption) however caused
// and on any theory of liability, whether in contract, strict liability,
// or tort (including negligence or otherwise) arising in any way out of
// the use of this software, even if advised of the possibility of such damage.
//
//M*/
#include "precomp.hpp"
#include "opencv2/imgproc.hpp"
using namespace std;
#define EPS 0.001F
namespace cv
{
namespace optflow
{
class DISOpticalFlowImpl : public DISOpticalFlow
{
public:
DISOpticalFlowImpl();
void calc(InputArray I0, InputArray I1, InputOutputArray flow);
void collectGarbage();
protected: // algorithm parameters
int finest_scale, coarsest_scale;
int patch_size;
int patch_stride;
int grad_descent_iter;
public: // getters and setters
int getFinestScale() const { return finest_scale; }
void setFinestScale(int val) { finest_scale = val; }
int getPatchSize() const { return patch_size; }
void setPatchSize(int val) { patch_size = val; }
int getPatchStride() const { return patch_stride; }
void setPatchStride(int val) { patch_stride = val; }
int getGradientDescentIterations() const { return grad_descent_iter; }
void setGradientDescentIterations(int val) { grad_descent_iter = val; }
private: // internal buffers
vector< Mat_<uchar> > I0s; // gaussian pyramid for the current frame
vector< Mat_<uchar> > I1s; // gaussian pyramid for the next frame
vector< Mat_<short> > I0xs; // gaussian pyramid for the x gradient of the current frame
vector< Mat_<short> > I0ys; // gaussian pyramid for the y gradient of the current frame
vector< Mat_<float> > Ux; // x component of the flow vectors
vector< Mat_<float> > Uy; // y component of the flow vectors
Mat_<Vec2f> U; // buffers for the merged flow
Mat_<float> Sx; // x component of the sparse flow vectors (before densification)
Mat_<float> Sy; // y component of the sparse flow vectors (before densification)
// structure tensor components and auxiliary buffers:
Mat_<float> I0xx_buf; // sum of squares of x gradient values
Mat_<float> I0yy_buf; // sum of squares of y gradient values
Mat_<float> I0xy_buf; // sum of x and y gradient products
Mat_<float> I0xx_buf_aux; // for computing sums using the summed area table
Mat_<float> I0yy_buf_aux;
Mat_<float> I0xy_buf_aux;
////////////////////////////////////////////////////////////
private: // private methods
void prepareBuffers(Mat &I0, Mat &I1);
void precomputeStructureTensor(Mat &dst_I0xx, Mat &dst_I0yy, Mat &dst_I0xy, Mat &I0x, Mat &I0y);
void patchInverseSearch(Mat &dst_Sx, Mat &dst_Sy, Mat &src_Ux, Mat &src_Uy, Mat &I0, Mat &I0x, Mat &I0y, Mat &I1);
void densification(Mat &dst_Ux, Mat &dst_Uy, Mat &src_Sx, Mat &src_Sy, Mat &I0, Mat &I1);
};
DISOpticalFlowImpl::DISOpticalFlowImpl()
{
finest_scale = 3;
patch_size = 9;
patch_stride = 4;
grad_descent_iter = 12;
}
void DISOpticalFlowImpl::prepareBuffers(Mat &I0, Mat &I1)
{
I0s.resize(coarsest_scale + 1);
I1s.resize(coarsest_scale + 1);
I0xs.resize(coarsest_scale + 1);
I0ys.resize(coarsest_scale + 1);
Ux.resize(coarsest_scale + 1);
Uy.resize(coarsest_scale + 1);
int fraction = 1;
int cur_rows = 0, cur_cols = 0;
for (int i = 0; i <= coarsest_scale; i++)
{
if (i == finest_scale)
{
cur_rows = I0.rows / fraction;
cur_cols = I0.cols / fraction;
I0s[i].create(cur_rows, cur_cols);
resize(I0, I0s[i], I0s[i].size(), 0.0, 0.0, INTER_AREA);
I1s[i].create(cur_rows, cur_cols);
resize(I1, I1s[i], I1s[i].size(), 0.0, 0.0, INTER_AREA);
Sx.create(cur_rows / patch_stride, cur_cols / patch_stride);
Sy.create(cur_rows / patch_stride, cur_cols / patch_stride);
I0xx_buf.create(cur_rows / patch_stride, cur_cols / patch_stride);
I0yy_buf.create(cur_rows / patch_stride, cur_cols / patch_stride);
I0xy_buf.create(cur_rows / patch_stride, cur_cols / patch_stride);
I0xx_buf_aux.create(cur_rows, cur_cols / patch_stride);
I0yy_buf_aux.create(cur_rows, cur_cols / patch_stride);
I0xy_buf_aux.create(cur_rows, cur_cols / patch_stride);
U.create(cur_rows, cur_cols);
}
else if (i > finest_scale)
{
cur_rows = I0s[i - 1].rows / 2;
cur_cols = I0s[i - 1].cols / 2;
I0s[i].create(cur_rows, cur_cols);
resize(I0s[i - 1], I0s[i], I0s[i].size(), 0.0, 0.0, INTER_AREA);
I1s[i].create(cur_rows, cur_cols);
resize(I1s[i - 1], I1s[i], I1s[i].size(), 0.0, 0.0, INTER_AREA);
}
fraction *= 2;
if (i >= finest_scale)
{
I0xs[i].create(cur_rows, cur_cols);
I0ys[i].create(cur_rows, cur_cols);
spatialGradient(I0s[i], I0xs[i], I0ys[i]);
Ux[i].create(cur_rows, cur_cols);
Uy[i].create(cur_rows, cur_cols);
}
}
}
void DISOpticalFlowImpl::precomputeStructureTensor(Mat &dst_I0xx, Mat &dst_I0yy, Mat &dst_I0xy, Mat &I0x, Mat &I0y)
{
short *I0x_ptr = I0x.ptr<short>();
short *I0y_ptr = I0y.ptr<short>();
float *I0xx_ptr = dst_I0xx.ptr<float>();
float *I0yy_ptr = dst_I0yy.ptr<float>();
float *I0xy_ptr = dst_I0xy.ptr<float>();
float *I0xx_aux_ptr = I0xx_buf_aux.ptr<float>();
float *I0yy_aux_ptr = I0yy_buf_aux.ptr<float>();
float *I0xy_aux_ptr = I0xy_buf_aux.ptr<float>();
int w = I0x.cols;
int h = I0x.rows;
int psz2 = patch_size / 2;
int psz = 2 * psz2;
// width of the sparse OF fields:
int ws = 1 + (w - patch_size) / patch_stride;
// separable box filter for computing patch sums on a sparse
// grid (determined by patch_stride)
for (int i = 0; i < h; i++)
{
float sum_xx = 0.0f, sum_yy = 0.0f, sum_xy = 0.0f;
short *x_row = I0x_ptr + i * w, *y_row = I0y_ptr + i * w;
for (int j = 0; j < patch_size; j++)
{
sum_xx += x_row[j] * x_row[j];
sum_yy += y_row[j] * y_row[j];
sum_xy += x_row[j] * y_row[j];
}
I0xx_aux_ptr[i * ws] = sum_xx;
I0yy_aux_ptr[i * ws] = sum_yy;
I0xy_aux_ptr[i * ws] = sum_xy;
int js = 1;
for (int j = patch_size; j < w; j++)
{
sum_xx += (x_row[j] * x_row[j] - x_row[j - patch_size] * x_row[j - patch_size]);
sum_yy += (y_row[j] * y_row[j] - y_row[j - patch_size] * y_row[j - patch_size]);
sum_xy += (x_row[j] * y_row[j] - x_row[j - patch_size] * y_row[j - patch_size]);
if ((j - psz) % patch_stride == 0)
{
I0xx_aux_ptr[i * ws + js] = sum_xx;
I0yy_aux_ptr[i * ws + js] = sum_yy;
I0xy_aux_ptr[i * ws + js] = sum_xy;
js++;
}
}
}
AutoBuffer<float> sum_xx_buf(ws), sum_yy_buf(ws), sum_xy_buf(ws);
float *sum_xx = (float *)sum_xx_buf;
float *sum_yy = (float *)sum_yy_buf;
float *sum_xy = (float *)sum_xy_buf;
for (int j = 0; j < ws; j++)
{
sum_xx[j] = 0.0f;
sum_yy[j] = 0.0f;
sum_xy[j] = 0.0f;
}
for (int i = 0; i < patch_size; i++)
for (int j = 0; j < ws; j++)
{
sum_xx[j] += I0xx_aux_ptr[i * ws + j];
sum_yy[j] += I0yy_aux_ptr[i * ws + j];
sum_xy[j] += I0xy_aux_ptr[i * ws + j];
}
for (int j = 0; j < ws; j++)
{
I0xx_ptr[j] = sum_xx[j];
I0yy_ptr[j] = sum_yy[j];
I0xy_ptr[j] = sum_xy[j];
}
int is = 1;
for (int i = patch_size; i < h; i++)
{
for (int j = 0; j < ws; j++)
{
sum_xx[j] += (I0xx_aux_ptr[i * ws + j] - I0xx_aux_ptr[(i - patch_size) * ws + j]);
sum_yy[j] += (I0yy_aux_ptr[i * ws + j] - I0yy_aux_ptr[(i - patch_size) * ws + j]);
sum_xy[j] += (I0xy_aux_ptr[i * ws + j] - I0xy_aux_ptr[(i - patch_size) * ws + j]);
}
if ((i - psz) % patch_stride == 0)
{
for (int j = 0; j < ws; j++)
{
I0xx_ptr[is * ws + j] = sum_xx[j];
I0yy_ptr[is * ws + j] = sum_yy[j];
I0xy_ptr[is * ws + j] = sum_xy[j];
}
is++;
}
}
}
void DISOpticalFlowImpl::patchInverseSearch(Mat &dst_Sx, Mat &dst_Sy, Mat &src_Ux, Mat &src_Uy, Mat &I0, Mat &I0x,
Mat &I0y, Mat &I1)
{
float *Ux_ptr = src_Ux.ptr<float>();
float *Uy_ptr = src_Uy.ptr<float>();
float *Sx_ptr = dst_Sx.ptr<float>();
float *Sy_ptr = dst_Sy.ptr<float>();
uchar *I0_ptr = I0.ptr<uchar>();
uchar *I1_ptr = I1.ptr<uchar>();
short *I0x_ptr = I0x.ptr<short>();
short *I0y_ptr = I0y.ptr<short>();
int w = I0.cols;
int h = I1.rows;
int psz2 = patch_size / 2;
// width and height of the sparse OF fields:
int ws = 1 + (w - patch_size) / patch_stride;
int hs = 1 + (h - patch_size) / patch_stride;
precomputeStructureTensor(I0xx_buf, I0yy_buf, I0xy_buf, I0x, I0y);
float *xx_ptr = I0xx_buf.ptr<float>();
float *yy_ptr = I0yy_buf.ptr<float>();
float *xy_ptr = I0xy_buf.ptr<float>();
// perform a fixed number of gradient descent iterations for each patch:
int i = psz2;
for (int is = 0; is < hs; is++)
{
int j = psz2;
for (int js = 0; js < ws; js++)
{
float cur_Ux = Ux_ptr[i * w + j];
float cur_Uy = Uy_ptr[i * w + j];
float detH = xx_ptr[is * ws + js] * yy_ptr[is * ws + js] - xy_ptr[is * ws + js] * xy_ptr[is * ws + js];
if (abs(detH) < EPS)
detH = EPS;
float invH11 = yy_ptr[is * ws + js] / detH;
float invH12 = -xy_ptr[is * ws + js] / detH;
float invH22 = xx_ptr[is * ws + js] / detH;
float prev_sum_diff = 100000000.0f;
for (int t = 0; t < grad_descent_iter; t++)
{
float dUx = 0, dUy = 0;
float diff = 0;
float sum_diff = 0.0f;
for (int pos_y = i - psz2; pos_y <= i + psz2; pos_y++)
for (int pos_x = j - psz2; pos_x <= j + psz2; pos_x++)
{
float pos_x_shifted = min(max(pos_x + cur_Ux, 0.0f), w - 1.0f - EPS);
float pos_y_shifted = min(max(pos_y + cur_Uy, 0.0f), h - 1.0f - EPS);
int pos_x_lower = (int)pos_x_shifted;
int pos_x_upper = pos_x_lower + 1;
int pos_y_lower = (int)pos_y_shifted;
int pos_y_upper = pos_y_lower + 1;
diff = (pos_x_shifted - pos_x_lower) * (pos_y_shifted - pos_y_lower) *
I1_ptr[pos_y_upper * w + pos_x_upper] +
(pos_x_upper - pos_x_shifted) * (pos_y_shifted - pos_y_lower) *
I1_ptr[pos_y_upper * w + pos_x_lower] +
(pos_x_shifted - pos_x_lower) * (pos_y_upper - pos_y_shifted) *
I1_ptr[pos_y_lower * w + pos_x_upper] +
(pos_x_upper - pos_x_shifted) * (pos_y_upper - pos_y_shifted) *
I1_ptr[pos_y_lower * w + pos_x_lower] -
I0_ptr[pos_y * w + pos_x];
sum_diff += diff * diff;
dUx += I0x_ptr[pos_y * w + pos_x] * diff;
dUy += I0y_ptr[pos_y * w + pos_x] * diff;
}
cur_Ux -= invH11 * dUx + invH12 * dUy;
cur_Uy -= invH12 * dUx + invH22 * dUy;
if (sum_diff > prev_sum_diff)
break;
prev_sum_diff = sum_diff;
}
if (norm(Vec2f(cur_Ux - Ux_ptr[i * w + j], cur_Uy - Uy_ptr[i * w + j])) <= patch_size)
{
Sx_ptr[is * ws + js] = cur_Ux;
Sy_ptr[is * ws + js] = cur_Uy;
}
else
{
Sx_ptr[is * ws + js] = Ux_ptr[i * w + j];
Sy_ptr[is * ws + js] = Uy_ptr[i * w + j];
}
j += patch_stride;
}
i += patch_stride;
}
}
void DISOpticalFlowImpl::densification(Mat &dst_Ux, Mat &dst_Uy, Mat &src_Sx, Mat &src_Sy, Mat &I0, Mat &I1)
{
float *Ux_ptr = dst_Ux.ptr<float>();
float *Uy_ptr = dst_Uy.ptr<float>();
float *Sx_ptr = src_Sx.ptr<float>();
float *Sy_ptr = src_Sy.ptr<float>();
uchar *I0_ptr = I0.ptr<uchar>();
uchar *I1_ptr = I1.ptr<uchar>();
int w = I0.cols;
int h = I0.rows;
int psz2 = patch_size / 2;
// width of the sparse OF fields:
int ws = 1 + (w - patch_size) / patch_stride;
int start_x;
int start_y;
start_y = psz2;
for (int i = 0; i < h; i++)
{
if (i - psz2 > start_y && start_y + patch_stride < h - psz2)
start_y += patch_stride;
start_x = psz2;
for (int j = 0; j < w; j++)
{
float coef, sum_coef = 0.0f;
float sum_Ux = 0.0f;
float sum_Uy = 0.0f;
if (j - psz2 > start_x && start_x + patch_stride < w - psz2)
start_x += patch_stride;
for (int pos_y = start_y; pos_y <= min(i + psz2, h - psz2 - 1); pos_y += patch_stride)
for (int pos_x = start_x; pos_x <= min(j + psz2, w - psz2 - 1); pos_x += patch_stride)
{
float diff;
int is = (pos_y - psz2) / patch_stride;
int js = (pos_x - psz2) / patch_stride;
float j_shifted = min(max(j + Sx_ptr[is * ws + js], 0.0f), w - 1.0f - EPS);
float i_shifted = min(max(i + Sy_ptr[is * ws + js], 0.0f), h - 1.0f - EPS);
int j_lower = (int)j_shifted;
int j_upper = j_lower + 1;
int i_lower = (int)i_shifted;
int i_upper = i_lower + 1;
diff = (j_shifted - j_lower) * (i_shifted - i_lower) * I1_ptr[i_upper * w + j_upper] +
(j_upper - j_shifted) * (i_shifted - i_lower) * I1_ptr[i_upper * w + j_lower] +
(j_shifted - j_lower) * (i_upper - i_shifted) * I1_ptr[i_lower * w + j_upper] +
(j_upper - j_shifted) * (i_upper - i_shifted) * I1_ptr[i_lower * w + j_lower] -
I0_ptr[i * w + j];
coef = 1 / max(1.0f, diff * diff);
sum_Ux += coef * Sx_ptr[is * ws + js];
sum_Uy += coef * Sy_ptr[is * ws + js];
sum_coef += coef;
}
Ux_ptr[i * w + j] = sum_Ux / sum_coef;
Uy_ptr[i * w + j] = sum_Uy / sum_coef;
}
}
}
void DISOpticalFlowImpl::calc(InputArray I0, InputArray I1, InputOutputArray flow)
{
CV_Assert(!I0.empty() && I0.depth() == CV_8U && I0.channels() == 1);
CV_Assert(!I1.empty() && I1.depth() == CV_8U && I1.channels() == 1);
CV_Assert(I0.sameSize(I1));
Mat I0Mat = I0.getMat();
Mat I1Mat = I1.getMat();
flow.create(I1Mat.size(), CV_32FC2);
Mat &flowMat = flow.getMatRef();
coarsest_scale = (int)(log((2 * I0Mat.cols) / (4.0 * patch_size))/log(2.0) + 0.5) - 1;
prepareBuffers(I0Mat, I1Mat);
Ux[coarsest_scale].setTo(0.0f);
Uy[coarsest_scale].setTo(0.0f);
for (int i = coarsest_scale; i >= finest_scale; i--)
{
patchInverseSearch(Sx, Sy, Ux[i], Uy[i], I0s[i], I0xs[i], I0ys[i], I1s[i]);
densification(Ux[i], Uy[i], Sx, Sy, I0s[i], I1s[i]);
// TODO: variational refinement step
if (i > finest_scale)
{
resize(Ux[i], Ux[i - 1], Ux[i - 1].size());
resize(Uy[i], Uy[i - 1], Uy[i - 1].size());
Ux[i - 1] *= 2;
Uy[i - 1] *= 2;
}
}
Mat uxy[] = {Ux[finest_scale], Uy[finest_scale]};
merge(uxy, 2, U);
resize(U, flowMat, flowMat.size());
flowMat *= pow(2, finest_scale);
}
void DISOpticalFlowImpl::collectGarbage()
{
I0s.clear();
I1s.clear();
I0xs.clear();
I0ys.clear();
Ux.clear();
Uy.clear();
U.release();
Sx.release();
Sy.release();
I0xx_buf.release();
I0yy_buf.release();
I0xy_buf.release();
I0xx_buf_aux.release();
I0yy_buf_aux.release();
I0xy_buf_aux.release();
}
Ptr<DISOpticalFlow> createOptFlow_DIS() { return makePtr<DISOpticalFlowImpl>(); }
}
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment