Browse Source

First version of slow_flow_release.

Joel Janai 1 year ago
parent
commit
48390ac151
48 changed files with 13886 additions and 1 deletions
  1. 131 0
      CMakeLists.txt
  2. 84 1
      README.md
  3. 698 0
      accumulate.cpp
  4. 695 0
      adaptiveFR.cpp
  5. 3 0
      adaptiveFR.dat
  6. 58 0
      cfgs/dense_tracking.cfg
  7. 60 0
      cfgs/slow_flow.cfg
  8. 41 0
      configuration.h
  9. 15 0
      configuration_epic.h
  10. 1961 0
      dense_tracking.cpp
  11. 48 0
      epic_flow_extended/Makefile
  12. 74 0
      epic_flow_extended/README
  13. 197 0
      epic_flow_extended/array_types.h
  14. 235 0
      epic_flow_extended/epic.cpp
  15. 31 0
      epic_flow_extended/epic.h
  16. 492 0
      epic_flow_extended/epic_aux.cpp
  17. 58 0
      epic_flow_extended/epic_aux.h
  18. BIN
      epic_flow_extended/epicflow-static
  19. 140 0
      epic_flow_extended/epicflow.cpp
  20. 791 0
      epic_flow_extended/image.c
  21. 114 0
      epic_flow_extended/image.h
  22. 403 0
      epic_flow_extended/io.c
  23. 28 0
      epic_flow_extended/io.h
  24. BIN
      epic_flow_extended/modelFinal.mat
  25. 399 0
      epic_flow_extended/solver.c
  26. 16 0
      epic_flow_extended/solver.h
  27. 144 0
      epic_flow_extended/variational.c
  28. 36 0
      epic_flow_extended/variational.h
  29. 302 0
      epic_flow_extended/variational_aux.c
  30. 35 0
      epic_flow_extended/variational_aux.h
  31. 926 0
      epic_flow_extended/variational_aux_mt.cpp
  32. 98 0
      epic_flow_extended/variational_aux_mt.h
  33. 843 0
      epic_flow_extended/variational_mt.cpp
  34. 73 0
      epic_flow_extended/variational_mt.h
  35. 17 0
      matlab/detect_edges.m
  36. 45 0
      penalty_functions/geman_mcclure.h
  37. 47 0
      penalty_functions/lorentzian.h
  38. 41 0
      penalty_functions/modified_l1_norm.h
  39. 22 0
      penalty_functions/penalty_function.h
  40. 31 0
      penalty_functions/quadratic_function.h
  41. 62 0
      penalty_functions/trunc_modified_l1_norm.h
  42. 1045 0
      slow_flow.cpp
  43. 519 0
      utils/hypothesis.cpp
  44. 269 0
      utils/hypothesis.h
  45. 823 0
      utils/parameter_list.cpp
  46. 145 0
      utils/parameter_list.h
  47. 1374 0
      utils/utils.cpp
  48. 217 0
      utils/utils.h

+ 131 - 0
CMakeLists.txt

@@ -0,0 +1,131 @@
1
+cmake_minimum_required(VERSION 2.8)
2
+
3
+project(slow_flow)
4
+
5
+set(CMAKE_C_COMPILER_INIT g++)
6
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11 -fPIC -Wall -g -O3 -msse4") 
7
+
8
+set(MIDDLEBURY_DEVKIT "[SPECIFY PATH TO MIDDLEBURY DEVKIT]")
9
+
10
+# opencv library
11
+find_package(OpenCV REQUIRED)
12
+
13
+# boost library
14
+find_package(Boost REQUIRED COMPONENTS system filesystem chrono regex)
15
+
16
+# eigen library
17
+find_package(Eigen3 REQUIRED )
18
+
19
+# PNG library
20
+find_package(PNG REQUIRED)
21
+
22
+# JPEG library
23
+find_package(JPEG REQUIRED)
24
+
25
+# Img library
26
+find_path(IMG_INCLUDE_DIRS Image.h ${MIDDLEBURY_DEVKIT}/cpp/imageLib)
27
+find_library(IMG_LIBRARY NAMES Img PATHS ${MIDDLEBURY_DEVKIT}/cpp/imageLib)
28
+add_library(Img STATIC IMPORTED)
29
+
30
+# LAPACK library
31
+find_package(LAPACK REQUIRED)
32
+find_package(BLAS REQUIRED)
33
+
34
+# GSL library
35
+FIND_PACKAGE(GSL REQUIRED)
36
+
37
+include_directories(${OpenCV_INCLUDE_DIRS})
38
+include_directories(${Eigen_INCLUDE_DIRS})
39
+include_directories(${Boost_INCLUDE_DIRS})
40
+include_directories(${PNG_INCLUDE_DIR})
41
+include_directories(${JPEG_INCLUDE_DIR})
42
+include_directories(${IMG_INCLUDE_DIRS})
43
+include_directories(${BLAS_INCLUDE_DIR})
44
+include_directories(${LAPACK_INCLUDE_DIR})
45
+include_directories(${GSLCBLAS_INCLUDE_DIR})
46
+include_directories(${GSL_INCLUDE_DIR})
47
+
48
+link_directories(${Boost_LIBRARY_DIR})
49
+
50
+# parallel computing
51
+find_package(OpenMP)
52
+if (OPENMP_FOUND)
53
+    message("OPENMP FOUND")
54
+    set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}")
55
+    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
56
+    set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_EXE_LINKER_FLAGS}")
57
+endif()
58
+
59
+add_executable(slow_flow slow_flow.cpp
60
+    epic_flow_extended/epic_aux.cpp
61
+    epic_flow_extended/epic.cpp
62
+    epic_flow_extended/image.c
63
+    epic_flow_extended/io.c
64
+    epic_flow_extended/solver.c
65
+    epic_flow_extended/variational_aux_mt.cpp
66
+    epic_flow_extended/variational_mt.cpp
67
+   	utils/utils.cpp
68
+    utils/parameter_list.cpp
69
+    ../libs/middlebury_devkit/cpp/colorcode.cpp
70
+    ../libs/middlebury_devkit/cpp/flowIO.cpp
71
+    ../libs/gco/GCoptimization.cpp
72
+    ../libs/gco/maxflow.cpp
73
+    ../libs/gco/graph.cpp
74
+    ../libs/gco/LinkedBlockList.cpp
75
+    ../libs/dmgunturk/dmha.c
76
+    ../libs/dmgunturk/dmbilinear.c
77
+    ../libs/dmgunturk/basic.c
78
+    ../libs/dmgunturk/imageio.c)
79
+
80
+add_executable(adaptiveFR adaptiveFR.cpp
81
+    epic_flow_extended/epic_aux.cpp
82
+    epic_flow_extended/epic.cpp
83
+    epic_flow_extended/image.c
84
+    epic_flow_extended/io.c
85
+    epic_flow_extended/solver.c
86
+    epic_flow_extended/variational_aux.c
87
+    epic_flow_extended/variational.c
88
+   	utils/utils.cpp
89
+    utils/parameter_list.cpp
90
+    ../libs/middlebury_devkit/cpp/colorcode.cpp
91
+    ../libs/middlebury_devkit/cpp/flowIO.cpp
92
+    ../libs/gco/GCoptimization.cpp
93
+    ../libs/gco/maxflow.cpp
94
+    ../libs/gco/graph.cpp
95
+    ../libs/gco/LinkedBlockList.cpp
96
+    ../libs/dmgunturk/dmha.c
97
+    ../libs/dmgunturk/dmbilinear.c
98
+    ../libs/dmgunturk/basic.c
99
+    ../libs/dmgunturk/imageio.c)
100
+
101
+add_executable(dense_tracking dense_tracking.cpp
102
+    epic_flow_extended/epic_aux.cpp
103
+    epic_flow_extended/epic.cpp
104
+    epic_flow_extended/image.c
105
+    epic_flow_extended/io.c
106
+    epic_flow_extended/solver.c
107
+    epic_flow_extended/variational_aux_mt.cpp
108
+    epic_flow_extended/variational_mt.cpp
109
+   	utils/hypothesis.cpp
110
+   	utils/utils.cpp
111
+    utils/parameter_list.cpp
112
+   	penalty_functions/penalty_function.h
113
+   	penalty_functions/lorentzian.h
114
+   	penalty_functions/modified_l1_norm.h
115
+   	penalty_functions/quadratic_function.h
116
+    ../libs/middlebury_devkit/cpp/colorcode.cpp
117
+    ../libs/middlebury_devkit/cpp/flowIO.cpp
118
+    ../libs/maxflow/maxflow.cpp
119
+    ../libs/maxflow/graph.cpp
120
+    ../libs/gco/GCoptimization.cpp
121
+    ../libs/gco/maxflow.cpp
122
+    ../libs/gco/graph.cpp
123
+    ../libs/gco/LinkedBlockList.cpp
124
+    ../libs/dmgunturk/dmha.c
125
+    ../libs/dmgunturk/dmbilinear.c
126
+    ../libs/dmgunturk/basic.c
127
+    ../libs/dmgunturk/imageio.c)
128
+
129
+target_link_libraries(slow_flow ${OpenCV_LIBS} ${Boost_LIBRARIES} ${IMG_LIBRARY} ${PNG_LIBRARY} ${JPEG_LIBRARY} ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES})
130
+target_link_libraries(adaptiveFR ${OpenCV_LIBS} ${Boost_LIBRARIES} ${IMG_LIBRARY} ${PNG_LIBRARY} ${JPEG_LIBRARY} ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES})
131
+target_link_libraries(dense_tracking ${OpenCV_LIBS} ${Boost_LIBRARIES} ${IMG_LIBRARY} ${PNG_LIBRARY} ${JPEG_LIBRARY} ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES} ${GSL_LIBRARY} ${GSLCBLAS_LIBRARY})

+ 84 - 1
README.md

@@ -1,2 +1,85 @@
1
-# slow_flow_release
1
+# Slow Flow: Generating Optical Flow Reference Data
2 2
 
3
+This code is based on the paper [Slow Flow: Exploiting High-Speed Cameras for
4
+Accurate and Diverse Optical Flow Reference Data](http://www.cvlibs.net/publications/Janai2017CVPR.pdf).
5
+
6
+The high speed flow code (slow_flow) is based on EpicFlow v1.00. We extended it to reason over multiple frames and occlusions. The extended code is included and for details on the orginal code we refer to the paper by Revaud et al. [EpicFlow: Edge-Preserving Interpolation of
7
+Correspondences for Optical Flow](https://hal.inria.fr/hal-01142656/document) and the project webpage 
8
+https://thoth.inrialpes.fr/src/epicflow/
9
+
10
+We provide [two teaser sequences]() to run our code. We are working on publishing the complete high speed datasets used in the project.
11
+
12
+## Compilation
13
+#### The following libraries are necessary to compile and use our code:
14
+
15
+	Eigen3, Boost, Atlas, Blas, Lapack, Flann, GSL, PNG and JPEG
16
+	sudo apt-get install libeigen3-dev libboost-all-dev libatlas-base-dev libblas-dev liblapack-dev libflann-dev libgsl-dev libpng-dev libjpeg-dev
17
+	 
18
+	Download and install opencv2.4
19
+	http://docs.opencv.org/2.4/doc/tutorials/introduction/linux_install/linux_install.html
20
+	
21
+	Download SED for edges
22
+	https://www.microsoft.com/en-us/download/details.aspx?id=52370&from=http%3A%2F%2Fresearch.microsoft.com%2Fen-us%2Fdownloads%2F389109f6-b4e8-404c-84bf-239f7cbf4e3d%2F
23
+
24
+	Download Piotr Dollar's toolbox
25
+	http://vision.ucsd.edu/~pdollar/toolbox/doc/index.html
26
+
27
+	Download Deep Matching 
28
+	(Optional: using course-to-fine by setting ‘deep_matching’ to 0 and ‘slow_flow_layers’  larger than 1)
29
+	http://lear.inrialpes.fr/people/revaud
30
+	
31
+	Download flow-code from Middlebury OF dataset and compile ImageLib
32
+	http://vision.middlebury.edu/flow/code/flow-code.zip
33
+	
34
+	Download gco-v3.0 library
35
+	http://vision.csd.uwo.ca/code/
36
+	
37
+	(Optional) Download Gunturk-Altunbasak-Mersereau Alternating Projections Image Demosaicking 
38
+	By setting ‘raw_demosaicing’ to 1 in the configuration file and uncommenting line 17 and 38 in configuration.h
39
+	http://www.ipol.im/pub/art/2011/g_gapd/
40
+
41
+#### The paths to libraries need to be specified in the following files:
42
+	configuration.h
43
+	configuration_epic.h
44
+	matlab/detect_edges.m
45
+	CMakeLists.txt 
46
+
47
+#### Run cmake and make to compile the code.
48
+
49
+## Run Pipeline ###
50
+#### 0. Run epic flow on low resolution to use adaptive frame rates for slow flow
51
+	./adaptiveFR -path [path] -folder [folder]
52
+
53
+	Examples for our teaser_sequences:
54
+	./adaptiveFR -path '[path to teaser]/sequence/' -folder 'sheeps' -raw
55
+	./adaptiveFR -path '[path to teaser]/sequence/' -folder 'ambush_2' -format 'out_%i_%03i.png' -start 491 -sintel
56
+
57
+#### 1. Run slow_flow for flow estimations of all consecutive high speed frames
58
+	./slow_flow [slow flow cfg file]
59
+
60
+	An example configuration file is provided in "cfgs".
61
+	
62
+	Optional: 
63
+		-jet		compute the flow for one specified high speed pair
64
+		-fr		using adaptive frame rate compute flow for 0: high frame rate, 1: low frame rate 
65
+		-resume		resume processing configuration file
66
+		-deep_settings	specify settings for deep matching
67
+
68
+#### 2. Run dense_tracking using the output of slow flow 
69
+	./dense_tracking [dense tracking cfg file]
70
+
71
+	Optional: 
72
+		-select		compute the flow for one specific final image pair
73
+		-resume		resume processing configuration file
74
+
75
+	An example configuration file is provided in "cfgs".
76
+
77
+## Citation
78
+
79
+If you use our code, please cite our paper:
80
+@INPROCEEDINGS{Janai2017CVPR,
81
+  author = {Joel Janai and Fatma Güney and Jonas Wulff and Michael Black and Andreas Geiger},
82
+  title = {Slow Flow: Exploiting High-Speed Cameras for Accurate and Diverse Optical Flow Reference Data},
83
+  booktitle = {Conference on Computer Vision and Pattern Recognition (CVPR)},
84
+  year = {2017}
85
+} 

+ 698 - 0
accumulate.cpp

@@ -0,0 +1,698 @@
1
+/*
2
+ * accumulate.cpp
3
+ *
4
+ *  Created on: Mar 7, 2016
5
+ *      Author: Janai
6
+ */
7
+
8
+#include <iostream>
9
+#include <fstream>
10
+#include <stdio.h>
11
+#include <string>
12
+#include <cmath>
13
+#include <omp.h>
14
+#include <random>
15
+
16
+#include <opencv2/core.hpp>
17
+#include <opencv2/highgui.hpp>
18
+#include <opencv2/imgproc.hpp>
19
+
20
+#include <boost/filesystem.hpp>
21
+
22
+#include "epic.h"
23
+#include "image.h"
24
+#include "io.h"
25
+#include "variational_mt.h"
26
+#include "utils/edge_detector.h"
27
+#include "utils/utils.h"
28
+
29
+#include "../multi_frame_flow/multi_frame_optical_flow.h"
30
+#include "../multi_frame_flow/utils/parameter_list.h"
31
+
32
+#include "../libs/middlebury_devkit/cpp/colorcode.h"
33
+#include "../libs/middlebury_devkit/cpp/flowIO.h"
34
+
35
+using namespace std;
36
+using namespace cv;
37
+
38
+void setDefaultVariational(ParameterList& params) {
39
+    // general
40
+    params.insert("verbose", "0", true);
41
+    params.insert("threads", "1", true);
42
+
43
+//    params.insert("format", "img_%04d.tif", true);
44
+
45
+    params.insert("stats", "0", true);
46
+    params.insert("S", "2", true);
47
+
48
+    params.insert("scale", "1.0f", true);
49
+
50
+    params.insert("dataterm", "1");
51
+    params.insert("smoothing", "1");
52
+
53
+    params.insert("epic_alpha", "1.0f");
54
+    params.insert("epic_gamma", "0.72f", true);
55
+    params.insert("epic_delta", "0.0f", true);
56
+
57
+    params.insert("layers", "1", true);
58
+    params.insert("p_scale", "0.9f", true);;
59
+
60
+    params.insert("niter_alter", "1", true);
61
+    params.insert("niter_outer", "5", true);
62
+    params.insert("thres_outer", "1e-5", true);
63
+    params.insert("niter_inner", "1", true);
64
+    params.insert("thres_inner", "1e-5", true);
65
+    params.insert("niter_solver", "30", true);
66
+    params.insert("sor_omega", "1.9f", true);
67
+
68
+    params.insert("robust_color", "1", true);
69
+    params.insert("robust_color_eps", "0.001", true);
70
+    params.insert("robust_color_truncation", "0.5", true);
71
+    params.insert("robust_reg", "1", true);
72
+    params.insert("robust_reg_eps", "0.001", true);
73
+    params.insert("robust_reg_truncation", "0.5", true);
74
+
75
+    params.insert("gradient_sigmoid", "5.0f", true);
76
+}
77
+
78
+int main(int argc, char **argv) {
79
+	/*
80
+	 *  getting config file location
81
+	 */
82
+#define isarg(key)  !strcmp(a,key)
83
+	string file;
84
+	string input_path = "";
85
+	string output_path = "";
86
+	if (argc > 1) {
87
+		file = string(argv[1]);
88
+
89
+		if(boost::filesystem::exists(file)) {
90
+			printf("using parameters %s\n", file.c_str());
91
+		} else {
92
+			cerr << "usage: ./accumulate [cfg-file] -input_path [path] -output_path [path] {optional:Frames}" << endl;
93
+			return -1;
94
+		}
95
+
96
+		int selected_jet = -1;
97
+		int current_arg = 1;
98
+		while(current_arg < argc ){
99
+			const char* a = argv[current_arg++];
100
+			if(a[0] != '-') {
101
+				continue;
102
+			}
103
+
104
+			if( isarg("-h") || isarg("-help") )
105
+				cerr << "usage: ./accumulate [cfg-file] -input_path [path] -output_path [path] {optional:Frames}" << endl;
106
+			else if( isarg("-input_path") )
107
+				input_path = string(argv[current_arg++]);
108
+			else if( isarg("-output_path") )
109
+				output_path = string(argv[current_arg++]);
110
+			else{
111
+				fprintf(stderr, "unknown argument %s", a);
112
+				cerr << "usage: ./accumulate [cfg-file] -input_path [path] -output_path [path] {optional:Frames}" << endl;
113
+				exit(1);
114
+			}
115
+		}
116
+	} else {
117
+		cerr << "usage: ./accumulate [cfg-file] -input_path [path] -output_path [path] {optional:Frames}" << endl;
118
+		return -1;
119
+	}
120
+
121
+	/*
122
+	 *  read in parameters and print config
123
+	 */
124
+	ParameterList params;
125
+
126
+	params.read(file);
127
+	string path = file.substr(0, file.find_last_of('/'));
128
+
129
+	bool sintel = params.parameter<bool>("sintel", "0");
130
+
131
+	// add input path and output path
132
+	params.setParameter("flow_file", input_path + params.parameter("flow_file"));
133
+	for(int i = 0; i < params.file_list.size(); i++)
134
+		params.file_list[i] = input_path + params.file_list[i];
135
+	for(int i = 0; i < params.file_gt_list.size(); i++)
136
+		params.file_gt_list[i] = input_path + params.file_gt_list[i];
137
+	params.output = path + "/";
138
+
139
+	/*
140
+	 * decompose sequence in subsets and adjust number of frames if necessary
141
+	 */
142
+	uint32_t Jets = params.Jets;
143
+
144
+	if (!params.exists("S") || params.parameter<int>("S") < 2)
145
+		params.setParameter<int>("S", 2);
146
+
147
+	int steps = params.parameter<int>("S") - 1;
148
+
149
+	if(Jets == 0) {
150
+		Jets = floor(params.F / steps);
151
+		params.Jets = Jets;
152
+	}
153
+
154
+
155
+//	int batch_size = floor((1.0f * params.parameter<int>("max_fps")) / (params.parameter<int>("ref_fps")));
156
+	int batch_size = floor((1.0f * params.parameter<int>("max_fps")) / (params.parameter<int>("jet_fps")));
157
+//	int batches = floor((1.0f * params.Jets) * batch_size);
158
+	Jets = floor(batch_size / steps);
159
+	int batches = params.parameter<int>("batches", "1");
160
+
161
+	int skip = floor((1.0f * params.parameter<int>("max_fps")) / (params.parameter<int>("jet_fps")));
162
+	if(params.exists("frames")) {
163
+		// decompose sequence in batches
164
+		skip = ceil(1.0f * batch_size / (params.parameter<int>("frames") - 1.0));			// with ceil we will estimate too far by fractions
165
+		Jets = floor(1.0f * Jets / skip);
166
+	}
167
+
168
+	// iterate over experiments
169
+	for(uint32_t exp = 0; exp < params.experiments(); exp++) {
170
+		ParameterList exp_params;
171
+
172
+		exp_params = ParameterList(params);
173
+		params.nextExp();
174
+
175
+		// get scaling factor
176
+		double scale = params.parameter<float>("scale", "1.0");
177
+
178
+		string flow_format = exp_params.parameter<string>("flow_format", "frame_%i");
179
+
180
+		// create folder for accumulated batches
181
+		stringstream acc_folder;
182
+		acc_folder << exp_params.output << "accumulation/";
183
+		boost::filesystem::create_directories(acc_folder.str());
184
+		boost::filesystem::create_directories(acc_folder.str() + "gt_occlusions/");          // accumulate folder
185
+
186
+		exp_params.insert("accumulation", acc_folder.str() , true);
187
+		exp_params.Jets = batches;
188
+		stringstream news;
189
+		news << batch_size*steps + 1;
190
+		exp_params.setParameter("S", news.str());
191
+
192
+		exp_params.print();
193
+
194
+		ofstream infos;
195
+		infos.open((acc_folder.str() + "config.cfg").c_str());
196
+		infos << "# SlowFlow variational estimation\n";
197
+		infos << exp_params;
198
+		infos.close();
199
+
200
+		uint32_t num_sources = params.file_list.size();
201
+
202
+		for (uint32_t source = 0; source < num_sources; source++) {
203
+			for(int b = 0; b < batches; b++) {
204
+				if(b > 0)
205
+					exp_params.sequence_start_list[source] += Jets * skip * steps;
206
+
207
+				// read in sequence, forward and backward flow, and segmentation
208
+				Mat* gt = NULL;  	 	 // plus reference frame
209
+				Mat* gt_occlusions = NULL;
210
+
211
+				Mat *forward_flow = new Mat[Jets];
212
+				Mat *backward_flow = new Mat[Jets];
213
+
214
+				/*
215
+				 * ################### read in ground truth ###################
216
+				 */
217
+				if(exp_params.file_gt.size() > 0) {
218
+					gt = new Mat[batch_size];
219
+					for (uint32_t f = 0; f < (batch_size); f++) {
220
+						char gtF[200];
221
+
222
+						if(!sintel)
223
+							sprintf(gtF, exp_params.file_gt_list[source].c_str(), exp_params.sequence_start_list[source] + f);
224
+						else {
225
+							int sintel_frame = exp_params.sequence_start_list[source] / 1000;
226
+							int hfr_frame = f + (exp_params.sequence_start_list[source] % 1000);
227
+
228
+							while(hfr_frame < 0) {
229
+								sintel_frame--;
230
+								hfr_frame = 42 + hfr_frame;
231
+							}
232
+							while(hfr_frame > 41) {
233
+								sintel_frame++;
234
+								hfr_frame = hfr_frame - 42;
235
+							}
236
+
237
+							sprintf(gtF, exp_params.file_gt_list[source].c_str(), sintel_frame, hfr_frame);
238
+						}
239
+
240
+						if (access(gtF, F_OK) == -1) {
241
+							continue;
242
+							cerr << "Error reading "<< gtF << "!" << endl;
243
+						}
244
+						cout << "Reading "<< gtF << "!" << endl;
245
+
246
+						gt[f] = readGTMiddlebury(gtF);
247
+
248
+						// rescale gt flow
249
+						resize(gt[f], gt[f], Size(0, 0), scale, scale, INTER_LINEAR);
250
+						gt[f] = scale * gt[f];
251
+					}
252
+				}
253
+
254
+				/*
255
+				 * ################### read in ground truth occlusions ###################
256
+				 */
257
+				if(exp_params.occlusions_list.size() > 0) {
258
+					gt_occlusions = new Mat[batch_size];
259
+					for (u_int32_t f = 0; f < (uint32_t) (batch_size); f++) {
260
+						char oocF[500];
261
+
262
+						if(!sintel)
263
+							sprintf(oocF, exp_params.occlusions_list[source].c_str(), exp_params.sequence_start_list[source] + f);
264
+						else {
265
+							int sintel_frame = exp_params.sequence_start_list[source] / 1000;
266
+							int hfr_frame = f + (exp_params.sequence_start_list[source] % 1000);
267
+
268
+							while(hfr_frame < 0) {
269
+								sintel_frame--;
270
+								hfr_frame = 42 + hfr_frame;
271
+							}
272
+							while(hfr_frame > 41) {
273
+								sintel_frame++;
274
+								hfr_frame = hfr_frame - 42;
275
+							}
276
+
277
+							sprintf(oocF, exp_params.occlusions_list[source].c_str(), sintel_frame, hfr_frame);
278
+						}
279
+
280
+						// read file only if already processed by first experiment
281
+						if(access( oocF, F_OK ) != -1) {
282
+							cout << "Reading "<< oocF << "!" << endl;
283
+
284
+							gt_occlusions[f] = imread(string(oocF));
285
+
286
+							resize(gt_occlusions[f], gt_occlusions[f], Size(0, 0), scale, scale, INTER_CUBIC);
287
+
288
+							memset(oocF, 0, 500);
289
+							sprintf(oocF, "%s/gt_occlusions/occ_%05i.png", acc_folder.str().c_str(), exp_params.sequence_start + f);
290
+							// read file only if already processed by first experiment
291
+							if(access( oocF, F_OK ) == -1) {
292
+								// write flow and flow image to file
293
+								imwrite(oocF, gt_occlusions[f]);
294
+							}
295
+
296
+							if(gt_occlusions[f].channels() > 1)
297
+								cvtColor(gt_occlusions[f], gt_occlusions[f], CV_BGR2GRAY);
298
+							gt_occlusions[f].convertTo(gt_occlusions[f], CV_8UC1);
299
+						} else {
300
+							cerr << "Error reading "<< oocF << "!" << endl;
301
+						}
302
+					}
303
+				}
304
+
305
+
306
+				char source_str[1024];
307
+				if(params.name_list.size() > 1)
308
+					sprintf(source_str, "%s%02i_", params.name_list[source].c_str(), params.id(source));
309
+				else
310
+					source_str[0] = '\0';
311
+
312
+				/*
313
+				 * ################### read in flow fields ###################
314
+				 */
315
+				for (uint32_t f = 0; f < Jets; f++) {
316
+					char fFlowF[200];
317
+					char bFlowF[200];
318
+
319
+	//				if (exp_params.exists("flow_file")) {
320
+	//					if(!sintel) {
321
+	//						sprintf(fFlowF, (exp_params.parameter("flow_file") + ".flo").c_str(), exp_params.sequence_start_list[source] + f * steps);
322
+	//						sprintf(bFlowF, (exp_params.parameter("flow_file") + "_back.flo").c_str(), exp_params.sequence_start_list[source] + f * steps + steps);
323
+	//					} else {
324
+	//						sprintf(fFlowF, (exp_params.parameter("flow_file") + ".flo").c_str(), exp_params.sequence_start_list[source] + f * steps, 0);
325
+	//						sprintf(bFlowF, (exp_params.parameter("flow_file") + "_back.flo").c_str(), exp_params.sequence_start_list[source] + f * steps + steps, 0);
326
+	//					}
327
+	//				} else {
328
+	//					sprintf(fFlowF, "%sframe_%i.flo", exp_params.output.c_str(), (exp_params.sequence_start_list[source] + f * steps));
329
+	//					sprintf(bFlowF, "%sframe_%i_back.flo", exp_params.output.c_str(), (exp_params.sequence_start_list[source] + f * steps + steps));
330
+	//				}
331
+					string hfr = "";
332
+					if (boost::filesystem::exists( exp_params.output + "hfr/" ))
333
+						hfr = "hfr/";
334
+
335
+					sprintf(fFlowF, ("%s%s" + flow_format + ".flo").c_str(), (exp_params.output + hfr).c_str(), source_str, (exp_params.sequence_start_list[source] + f * steps * skip));
336
+					if (exp_params.parameter<bool>("wrong_backward", "0"))
337
+						sprintf(bFlowF, ("%s%s" + flow_format + "_back.flo").c_str(), (exp_params.output + hfr).c_str(), source_str,
338
+								(exp_params.sequence_start_list[source] + f * steps * skip));
339
+					else
340
+						sprintf(bFlowF, ("%s%s" + flow_format + "_back.flo").c_str(), (exp_params.output + hfr).c_str(), source_str,
341
+								(exp_params.sequence_start_list[source] + f * steps * skip + steps * skip));
342
+
343
+					if (!boost::filesystem::exists(fFlowF)) {
344
+						cerr << fFlowF << " does not exist!" << endl;
345
+						return -1;
346
+					}
347
+					if (!boost::filesystem::exists(bFlowF)) {
348
+						cerr << bFlowF << " does not exist!" << endl;
349
+						return -1;
350
+					}
351
+
352
+					cout << "Reading " << fFlowF <<  endl;
353
+					cout << "Reading " << bFlowF <<  endl;
354
+
355
+					forward_flow[f] = readGTMiddlebury(fFlowF);
356
+					backward_flow[f] = readGTMiddlebury(bFlowF);
357
+
358
+					// crop images
359
+					if(exp_params.center.x > 0) {
360
+						// NOTE: ROWRANGE IS INDUCING ERRORS IN ACCUMULATION!!!!!!
361
+						forward_flow[f] = crop(forward_flow[f], exp_params.center, exp_params.extent);
362
+						backward_flow[f] = crop(backward_flow[f], exp_params.center, exp_params.extent);
363
+					}
364
+				}
365
+				Mat lfr_occlusions = Mat::zeros(forward_flow[0].rows, forward_flow[0].cols, CV_8UC1);
366
+				if(gt != NULL && gt_occlusions != NULL) lfr_occlusions = fuseOcclusions(gt, gt_occlusions, 0, Jets);
367
+
368
+				// accumulate ground truth flow
369
+				Mat fconf, bconf;
370
+				Mat* acc_flow = new Mat[Jets];
371
+//				accumulateFlow(acc_flow,  &forward_flow[b * Jets], lfr_occlusions, Jets);
372
+				accumulateConsistentBatches(acc_flow, forward_flow, backward_flow, fconf, bconf, NULL, Jets, params.parameter<double>("acc_consistency_threshold_slope"), 0, params.parameter<bool>("acc_discard_inconsistent", "1"), 0);
373
+
374
+				/*
375
+				 *  ########################################### Refine flow ###########################################
376
+				 */
377
+
378
+//				if(exp_params.parameter<bool>("refine", "1")) {
379
+//					int width = acc_flow->cols;
380
+//					int height = acc_flow->rows;
381
+//
382
+//					ParameterList refine_params(exp_params);
383
+//					setDefaultVariational(refine_params);
384
+//					refine_params.F = 2;
385
+//					refine_params.Jets = 1;
386
+//					refine_params.insert("S", "2", true);
387
+//					refine_params.verbose = "1000";
388
+//					refine_params.insert("layers", "1", true);
389
+//					refine_params.insert("method", "forward", true);
390
+//					refine_params.insert("niter_alter", "1", true);
391
+//					refine_params.insert("niter_outer", "5", true);
392
+//					refine_params.insert("occlusion_reasoning", "0", true);
393
+//					refine_params.insert("framerate_reasoning", "0", true);
394
+//					refine_params.insert("rho_0", "1", true);
395
+//					refine_params.insert("omega_0", "0", true);
396
+//
397
+//					color_image_t** im = new color_image_t*[3];
398
+//					im[0] = color_image_new(width, height);
399
+//					im[1] = color_image_new(width, height);
400
+//					im[2] = color_image_new(width, height);
401
+//
402
+
403
+//					/*
404
+//					 * ################### read in image sequence ###################
405
+//					 */
406
+//					for (uint32_t f = 0; f < 2; f++) {
407
+//						char img_file[1024];
408
+//						if (!sintel) {
409
+//							sprintf(img_file, (sequence_path + format).c_str(), thread_params.sequence_start_list[source] + f * Jets * steps * skip);
410
+//						} else {
411
+//							int sintel_frame = thread_params.sequence_start_list[source] / 1000;
412
+//							int hfr_frame = f * steps * skip + (thread_params.sequence_start_list[source] % 1000);
413
+//
414
+//							while (hfr_frame < 0) {
415
+//								sintel_frame--;
416
+//								hfr_frame = 42 + hfr_frame;
417
+//							}
418
+//							while (hfr_frame > 41) {
419
+//								sintel_frame++;
420
+//								hfr_frame = hfr_frame - 42;
421
+//							}
422
+//
423
+//							sprintf(img_file, (sequence_path + format).c_str(), sintel_frame, hfr_frame);
424
+//						}
425
+//						cout << "Reading " << img_file << "..." << endl;
426
+//
427
+//						sequence[f] = imread(string(img_file), CV_LOAD_IMAGE_UNCHANGED);         // load images
428
+//
429
+//						float norm = 1;
430
+//						if (sequence[f].type() == 2 || sequence[f].type() == 18)
431
+//							norm = 1.0f / 255;		// for 16 bit images
432
+//
433
+//						// convert to floating point
434
+//						sequence[f].convertTo(sequence[f], CV_32FC(sequence[f].channels()));
435
+//
436
+//						/*
437
+//						 * DEMOSAICING
438
+//						 */
439
+//						if (thread_params.exists("raw") && thread_params.parameter<bool>("raw")) {
440
+//							Mat tmp = sequence[f].clone();
441
+//							color_image_t* tmp_in = color_image_new(sequence[f].cols, sequence[f].rows);
442
+//							color_image_t* tmp_out = color_image_new(sequence[f].cols, sequence[f].rows);
443
+//
444
+//							switch (thread_params.parameter<int>("raw_demosaicing", "0")) {
445
+//								case 0: // use bilinear demosaicing
446
+//									sequence[f] = Mat::zeros(tmp.rows, tmp.cols, CV_32FC3);
447
+//									//						bayer2rgb(tmp, sequence[f], green_start, blue_start); 	// red green
448
+//									bayer2rgbGR(tmp, sequence[f], red_loc[0], red_loc[1]); // red green
449
+//									break;
450
+//
451
+//								case 1:	// use hamilton adams demosaicing
452
+//									mat2colorImg<float>(sequence[f], tmp_in);
453
+//
454
+//									HamiltonAdamsDemosaic(tmp_out->c1, tmp_in->c1, tmp_in->width, tmp_in->height, red_loc[0], red_loc[1]); // Hamilton-Adams implemented by Pascal Getreuer
455
+//
456
+//									sequence[f] = Mat::zeros(sequence[f].rows, sequence[f].cols, CV_32FC3);
457
+//									colorImg2colorMat<Vec3f>(tmp_out, sequence[f]);
458
+//									break;
459
+//
460
+//								case 2: // use opencv demosaicing
461
+//									tmp.convertTo(tmp, CV_8UC1);
462
+//									sequence[f] = Mat::zeros(tmp.rows, tmp.cols, CV_8UC3);
463
+//
464
+//									int code = CV_BayerBG2RGB;
465
+//									if (red_loc[1] == 0) // y
466
+//										if (red_loc[0] == 0) // x
467
+//											code = CV_BayerBG2RGB;
468
+//										else
469
+//											code = CV_BayerGB2RGB;
470
+//									else if (red_loc[0] == 0) // x
471
+//										code = CV_BayerGR2RGB;
472
+//									else
473
+//										code = CV_BayerRG2RGB;
474
+//
475
+//									cv::cvtColor(tmp, sequence[f], code); // components from second row, second column !!!!!!!!!!!!!!!!!
476
+//									sequence[f].convertTo(sequence[f], CV_32FC(sequence[f].channels()));
477
+//									break;
478
+//							}
479
+//
480
+//							color_image_delete(tmp_in);
481
+//							color_image_delete(tmp_out);
482
+//						} else {
483
+//							// covert to RGB
484
+//							cv::cvtColor(sequence[f], sequence[f], CV_BGR2RGB);
485
+//						}
486
+//
487
+//						/*
488
+//						 * COLOR CORRECTION
489
+//						 */
490
+//						vector<Mat> channels;
491
+//						split(sequence[f], channels);
492
+//						float mid_val = (0.5 * 255.0 / norm);
493
+//						channels[0] = contrast * (channels[0] + red_balance + brightness - mid_val) + mid_val;
494
+//						channels[1] = contrast * (channels[1] + green_balance + brightness - mid_val) + mid_val;
495
+//						channels[2] = contrast * (channels[2] + blue_balance + brightness - mid_val) + mid_val;
496
+//						merge(channels, sequence[f]);
497
+//
498
+//						if (thread_params.parameter<bool>("grayscale", "0"))
499
+//							cvtColor(sequence[f], sequence[f], CV_RGB2GRAY);
500
+//
501
+//						// use only a part of the images
502
+//						if (thread_params.extent.x > 0 || thread_params.extent.y > 0) {
503
+//							sequence[f] = sequence[f].rowRange(Range(thread_params.center.y - thread_params.extent.y / 2, thread_params.center.y + thread_params.extent.y / 2));
504
+//							sequence[f] = sequence[f].colRange(Range(thread_params.center.x - thread_params.extent.x / 2, thread_params.center.x + thread_params.extent.x / 2));
505
+//						}
506
+//
507
+//						// rescale image with gaussian blur to avoid anti-aliasing
508
+//						double img_scale = thread_params.parameter<double>("scale", "1.0");
509
+//						if (img_scale != 1) {
510
+//							GaussianBlur(sequence[f], sequence[f], Size(0, 0), 1 / sqrt(2 * img_scale), 1 / sqrt(2 * img_scale), BORDER_REPLICATE);
511
+//							resize(sequence[f], sequence[f], Size(0, 0), img_scale, img_scale, INTER_LINEAR);
512
+//						}
513
+//
514
+//						// print to file
515
+//						char file[1024];
516
+//						sprintf(file, (acc_folder.str() + "sequence/frame_%s%i.png").c_str(), source_str, thread_params.sequence_start_list[source] - steps * skip + f * skip);
517
+//
518
+//						Mat output_img;
519
+//						if (thread_params.parameter<bool>("16bit", "0")) {
520
+//							sequence[f].convertTo(output_img, CV_16UC(sequence[f].channels()));
521
+//						} else {
522
+//							sequence[f].convertTo(output_img, CV_8UC(sequence[f].channels()), norm);
523
+//						}
524
+//
525
+//						if (thread_params.parameter<bool>("grayscale", "0"))
526
+//							cv::cvtColor(output_img, output_img, CV_GRAY2BGR);	// OpenCV uses BGR
527
+//						else
528
+//							cv::cvtColor(output_img, output_img, CV_RGB2BGR);	// OpenCV uses BGR
529
+//
530
+//						if (thread_params.verbosity(WRITE_FILES)) {
531
+//							imwrite(file, output_img, compression_params);
532
+//						}
533
+//
534
+//						data[f] = color_image_new(sequence[f].cols, sequence[f].rows);
535
+//						if (thread_params.parameter<bool>("grayscale", "0"))
536
+//							mat2colorImg<float>(sequence[f], data[f]);
537
+//						else
538
+//							colorMat2colorImg<Vec3f>(sequence[f], data[f]);
539
+//
540
+//						if (thread_params.verbosity(VER_IN_GT) && (int) f < show_frames) {
541
+//							if (source == show_frames_source) {
542
+//								stringstream ftitle;
543
+//								ftitle << "Frame " << source;
544
+//								namedWindow(ftitle.str());               // Create a window for display.
545
+//								imshow(ftitle.str(), output_img);
546
+//
547
+//								// set the callback function for pixel selection
548
+//								int* cbparams = new int[2];
549
+//								cbparams[0] = f;
550
+//								cbparams[1] = source;
551
+//								setMouseCallback(ftitle.str(), CallBackFunc, cbparams);
552
+//
553
+//								waitKey(0);
554
+//							}
555
+//						}
556
+//					}
557
+//
558
+//					// normalize data terms
559
+//					normalize(data, Jets + 1, thread_params);
560
+//
561
+//					// rescale image with gaussian blur to avoid anti-aliasing
562
+//					double rescale = 1.0f * width / width;
563
+//					Mat tmp_1 = sequence[0].clone(),
564
+//						tmp_2 = sequence[Jets].clone();
565
+//					if (rescale != 1) {
566
+//						GaussianBlur(tmp_1, tmp_1, Size(0, 0), 1 / sqrt(2 * rescale), 1 / sqrt(2 * rescale), BORDER_REPLICATE);
567
+//						GaussianBlur(tmp_2, tmp_2, Size(0, 0), 1 / sqrt(2 * rescale), 1 / sqrt(2 * rescale), BORDER_REPLICATE);
568
+//					}
569
+//					resize(tmp_1, tmp_1, Size(width, height), 0, 0, INTER_LINEAR);
570
+//					resize(tmp_2, tmp_2, Size(width, height), 0, 0, INTER_LINEAR);
571
+//					if (thread_params.parameter<bool>("grayscale", "0")) {
572
+//						mat2colorImg<float>(tmp_2, im[0]);// assumes backward frame for symmetric option
573
+//						mat2colorImg<float>(tmp_1, im[1]);
574
+//						mat2colorImg<float>(tmp_2, im[2]);
575
+//					} else {
576
+//						colorMat2colorImg<Vec3f>(tmp_2, im[0]);// assumes backward frame for symmetric option
577
+//						colorMat2colorImg<Vec3f>(tmp_1, im[1]);
578
+//						colorMat2colorImg<Vec3f>(tmp_2, im[2]);
579
+//					}
580
+//
581
+//					normalize(im, 3, refine_params);
582
+//
583
+//					image_t *wx = image_new(im[0]->width, im[0]->height), *wy = image_new(im[0]->width, im[0]->height);
584
+//					vector<Mat> vu;
585
+//					split(flow[0][Jets - 1], vu);
586
+//					mat2img<double>(vu[1], wx);
587
+//					mat2img<double>(vu[0], wy);
588
+//
589
+////						bool change = true;
590
+////						while(change) {
591
+////							for(int i=0 ; i<wx->height ; i++){
592
+////								for( int j=0 ; j<wx->width ; j++){
593
+////									change = false;
594
+////
595
+////									if(wx->data[i*wx->stride+j] > UNKNOWN_FLOW_THRESH) {
596
+////										wx->data[i*wx->stride+j] = 0;
597
+////										int count = 0;
598
+////										if(wx->data[(i+1)*wx->stride+j] < UNKNOWN_FLOW_THRESH) {
599
+////											wx->data[i*wx->stride+j] += wx->data[(i+1)*wx->stride+j];
600
+////											count++;
601
+////										}
602
+////										if(wx->data[(i-1)*wx->stride+j] < UNKNOWN_FLOW_THRESH) {
603
+////											wx->data[i*wx->stride+j] += wx->data[(i-1)*wx->stride+j];
604
+////											count++;
605
+////										}
606
+////										if(wx->data[i*wx->stride+(j+1)] < UNKNOWN_FLOW_THRESH) {
607
+////											wx->data[i*wx->stride+j] += wx->data[i*wx->stride+(j+1)];
608
+////											count++;
609
+////										}
610
+////										if(wx->data[i*wx->stride+(j-1)] < UNKNOWN_FLOW_THRESH) {
611
+////											wx->data[i*wx->stride+j] += wx->data[i*wx->stride+(j-1)];
612
+////											count++;
613
+////										}
614
+////
615
+////										wx->data[i*wx->stride+j] = wx->data[i*wx->stride+j] / count;
616
+////										change = true;
617
+////									}
618
+////
619
+////									if(wy->data[i*wy->stride+j] > UNKNOWN_FLOW_THRESH) {
620
+////										wy->data[i*wy->stride+j] = 0;
621
+////										int count = 0;
622
+////										if(wy->data[(i+1)*wy->stride+j] < UNKNOWN_FLOW_THRESH) {
623
+////											wy->data[i*wy->stride+j] += wy->data[(i+1)*wy->stride+j];
624
+////											count++;
625
+////										}
626
+////										if(wy->data[(i-1)*wy->stride+j] < UNKNOWN_FLOW_THRESH) {
627
+////											wy->data[i*wy->stride+j] += wy->data[(i-1)*wy->stride+j];
628
+////											count++;
629
+////										}
630
+////										if(wy->data[i*wy->stride+(j+1)] < UNKNOWN_FLOW_THRESH) {
631
+////											wy->data[i*wy->stride+j] += wy->data[i*wy->stride+(j+1)];
632
+////											count++;
633
+////										}
634
+////										if(wy->data[i*wy->stride+(j-1)] < UNKNOWN_FLOW_THRESH) {
635
+////											wy->data[i*wy->stride+j] += wy->data[i*wy->stride+(j-1)];
636
+////											count++;
637
+////										}
638
+////
639
+////										wy->data[i*wy->stride+j] = wy->data[i*wy->stride+j] / count;
640
+////										change = true;
641
+////									}
642
+////								}
643
+////							}
644
+////						}
645
+//
646
+//					// energy minimization
647
+//				//	time(&t_start);
648
+//					Variational_MT minimzer_f;
649
+//					minimzer_f.variational(wx, wy, im, refine_params);
650
+//				//	time(&t_end);
651
+//
652
+//					img2mat<double>(wx, vu[1]);
653
+//					img2mat<double>(wy, vu[0]);
654
+//					merge(vu, flow[0][Jets - 1]);
655
+//
656
+//					color_image_delete(im[0]);
657
+//					color_image_delete(im[1]);
658
+//					color_image_delete(im[2]);
659
+//					delete[] im;
660
+//					image_delete(wx);
661
+//					image_delete(wy);
662
+//				}
663
+
664
+				// write final estimation
665
+//				stringstream flowF;
666
+//				flowF << acc_folder.str() << "/" << exp_params.name_list[source].c_str() << source  << "frame_" << exp_params.sequence_start_list[source] + steps * batch_size * b;
667
+				char flowF[200];
668
+				sprintf(flowF, ("%s/%s" + flow_format).c_str(), acc_folder.str().c_str(), source_str, (exp_params.sequence_start_list[source]));
669
+
670
+
671
+				Mat flow_img = flowColorImg(acc_flow[Jets-1], exp_params.verbosity(VER_CMD));
672
+
673
+				if (!flow_img.data) {                 // Check for invalid input
674
+					cout << "No flow for frame " << exp_params.sequence_start_list[source] << std::endl;
675
+					continue;
676
+				}
677
+
678
+				imwrite((string(flowF) + "_vis.png"), flow_img);
679
+
680
+				writeFlowMiddlebury(acc_flow[Jets-1], (string(flowF) + ".flo"));
681
+
682
+				//----------------------------------
683
+				// ########################################### final clean up ###########################################
684
+				//----------------------------------
685
+
686
+
687
+				if(gt != NULL) 				delete[] gt;
688
+				if(gt_occlusions != NULL) 	delete[] gt_occlusions;
689
+				delete[] forward_flow;
690
+				delete[] backward_flow;
691
+			}
692
+		}
693
+	}
694
+
695
+	return 0;
696
+}
697
+
698
+

+ 695 - 0
adaptiveFR.cpp

@@ -0,0 +1,695 @@
1
+/*
2
+ * adaptiveFR.cpp
3
+ *
4
+ *  Created on: Sep 15, 2016
5
+ *      Author: jjanai
6
+ */
7
+#include <fstream>
8
+#include <stdlib.h>
9
+#include <string.h>
10
+#include <string>
11
+#include <cmath>
12
+#include <omp.h>
13
+#include <unistd.h>
14
+#include <stdio.h>
15
+
16
+#include <boost/filesystem.hpp>
17
+#include <boost/regex.hpp>
18
+
19
+#include <opencv2/core.hpp>
20
+#include <opencv2/highgui.hpp>
21
+#include <opencv2/imgproc.hpp>
22
+
23
+#include "epic_flow_extended/image.h"
24
+#include "epic_flow_extended/io.h"
25
+#include "epic_flow_extended/epic.h"
26
+#include "epic_flow_extended/variational.h"
27
+#include "utils/utils.h"
28
+#include "utils/parameter_list.h"
29
+#include "configuration.h"
30
+
31
+using namespace std;
32
+using namespace cv;
33
+namespace fs = boost::filesystem;
34
+
35
+enum COMPARISON {GROUNDTRUTH = 0, WARPING = 1};
36
+
37
+/* show usage information */
38
+void usage(){
39
+    printf("usage:\n");
40
+    printf("    ./adaptiveFR -path [path] -folder [specific folder | file with list] -format [file format] -scale [default 0.25] -skip [target frame (2)] -sample [number of estimation (10)] -step [frames between estimation (10)] -start [first frame (0)] -quantil -raw -overwrite -sintel -subframes -threads\n");
41
+    printf("\n");
42
+}
43
+
44
+void setDefault(ParameterList& params) {
45
+	// general
46
+	    params.insert("verbose", "0", true);
47
+	    params.insert("threads", "1", true);
48
+
49
+	    params.insert("scale", "1.0f", true);
50
+
51
+	    params.insert("slow_flow_S", "2", true);
52
+
53
+	    // energy function
54
+	    params.insert("slow_flow_alpha", "4.0f");
55
+	    params.insert("slow_flow_gamma", "6.0f", true);
56
+	    params.insert("slow_flow_delta", "1.0f", true);
57
+
58
+	    // image pyramid
59
+	    params.insert("slow_flow_layers", "1", true);
60
+	    params.insert("slow_flow_p_scale", "0.9f", true);
61
+
62
+	    // optimization
63
+	    params.insert("slow_flow_niter_alter", "10", true);
64
+	    params.insert("slow_flow_niter_outer", "10", true);
65
+	    params.insert("slow_flow_thres_outer", "1e-5", true);
66
+	    params.insert("slow_flow_niter_inner", "1", true);
67
+	    params.insert("slow_flow_thres_inner", "1e-5", true);
68
+	    params.insert("slow_flow_niter_solver", "30", true);
69
+	    params.insert("slow_flow_sor_omega", "1.9f", true);
70
+
71
+	    // occlusion reasoning
72
+	    params.insert("slow_flow_occlusion_reasoning", "1", true);
73
+	    params.insert("slow_flow_occlusion_penalty", "0.1", true);
74
+	    params.insert("slow_flow_occlusion_alpha", "0.1", true);
75
+	    params.insert("slow_flow_output_occlusions", "1", true);
76
+
77
+	    // regularization
78
+	    params.insert("slow_flow_robust_color", "1", true);
79
+	    params.insert("slow_flow_robust_color_eps", "0.001", true);
80
+	    params.insert("slow_flow_robust_color_truncation", "0.5", true);
81
+	    params.insert("slow_flow_robust_reg", "1", true);
82
+	    params.insert("slow_flow_robust_reg_eps", "0.001", true);
83
+	    params.insert("slow_flow_robust_reg_truncation", "0.5", true);
84
+}
85
+
86
+inline bool insideImg(double x, double y, int width, int height) {
87
+	return (y >= 0 && y < height && x >= 0 && x < width);
88
+}
89
+
90
+int main(int argc, char **argv){
91
+    if( argc < 2){
92
+        if(argc>1) fprintf(stderr,"Error, not enough arguments\n");
93
+        usage();
94
+        exit(1);
95
+    }
96
+
97
+    // read optional arguments
98
+    string format = "%07i.tif";
99
+    uint32_t start = 0;
100
+
101
+	#define isarg(key)  !strcmp(a,key)
102
+	string path = "";
103
+	string folder = "";
104
+	string append = "";
105
+	bool overwrite = false;
106
+
107
+	int samples = 40;				// number of optical flow estimates
108
+	int sample_step = 10;			// step size
109
+	uint32_t all_frames = 2;		// frames for optical flow estimation (classical two frame formulation)
110
+	int skip = 2;
111
+
112
+	float q = 0.90f;				// quantil
113
+
114
+	double scale = 0.25;
115
+	double dm_scale = 1.0f;
116
+
117
+	bool sintel = false;			// specific file names (we would like to be able to distinguish frame number from 24 fps and 1008 fps)
118
+	bool subframes = false;			// are subframes specified
119
+
120
+	bool raw = false;
121
+
122
+	int threads = 1;
123
+
124
+	int current_arg = 0;
125
+	while(current_arg < argc ){
126
+		const char* a = argv[current_arg++];
127
+		if(a[0] != '-') {
128
+			continue;
129
+		}
130
+
131
+
132
+		if( isarg("-h") || isarg("-help") )
133
+			usage();
134
+		else if( isarg("-path") )
135
+			path = string(argv[current_arg++]);
136
+		else if( isarg("-folder") )
137
+			folder = string(argv[current_arg++]);
138
+		else if( isarg("-threads") )
139
+			threads = atoi(argv[current_arg++]);
140
+		else if( isarg("-append") )
141
+			append = string(argv[current_arg++]);
142
+		else if( isarg("-scale") )
143
+			scale = atof(argv[current_arg++]);
144
+		else if( isarg("-skip") )
145
+			skip = max(1, atoi(argv[current_arg++]));
146
+		else if( isarg("-samples") )
147
+			samples = atoi(argv[current_arg++]);
148
+		else if( isarg("-step") )
149
+			sample_step = atoi(argv[current_arg++]);
150
+		else if( isarg("-start") )
151
+			start = atoi(argv[current_arg++]);
152
+		else if( isarg("-quantil") )
153
+			q = atof(argv[current_arg++]);
154
+		else if( isarg("-overwrite") )
155
+			overwrite = true;
156
+		else if( isarg("-sintel") )
157
+			sintel = true;
158
+		else if( isarg("-raw") )
159
+			raw = true;
160
+		else if( isarg("-subframes") )
161
+			subframes = true;
162
+		else if( isarg("-format") )
163
+			format = string(argv[current_arg++]);
164
+		else {
165
+			fprintf(stderr, "unknown argument %s", a);
166
+			usage();
167
+			exit(1);
168
+		}
169
+	}
170
+
171
+	vector<string> folders;
172
+	if(folder.empty()) {
173
+		fs::path apk_path(path + "/");
174
+		boost::filesystem::directory_iterator dir(path + "/"), it, end;
175
+
176
+		for(it = dir; it != end; it++)
177
+		{
178
+			const boost::filesystem::path& p = *it;
179
+			folder = p.filename().string();
180
+			if(boost::filesystem::is_directory(p) &&
181
+				folder.compare("$RECYCLE.BIN") != 0 && folder.compare("preview") != 0 && folder.compare("Rallye") != 0 && folder.compare("System Volume Information") != 0 &&
182
+				folder.compare("WDApps") != 0 && folder.c_str()[0] != '.') {
183
+
184
+				folders.push_back(p.filename().string());
185
+			}
186
+
187
+		}
188
+	} else {
189
+		if(boost::filesystem::is_directory(path + "/" + folder + "/")) {
190
+			folders.push_back(folder);
191
+		} else {
192
+			ifstream folder_input;
193
+			string ifilename(folder);
194
+			if (std::strcmp(ifilename.c_str(), "-") != 0) {
195
+				folder_input.open(ifilename.c_str());
196
+
197
+				if (!folder_input.is_open()) {
198
+					std::cerr << ifilename << ": " << "no such file or directory" << "\n";
199
+					return EXIT_FAILURE;
200
+				}
201
+
202
+				string line;
203
+
204
+				// parse header
205
+				while(getline(folder_input, line)) {
206
+					if(boost::filesystem::is_directory(path + "/" + line + "/")) {
207
+						folders.push_back(line);
208
+					} else
209
+						std::cerr << path + "/" + line + "/" << ": " << "no such directory" << "\n";
210
+				}
211
+			}
212
+			folder_input.close();
213
+		}
214
+	}
215
+
216
+	sort(folders.begin(), folders.end());
217
+
218
+	if(sintel && !subframes)
219
+		start = start * 1000;
220
+
221
+	stringstream overview;
222
+	#pragma omp parallel for num_threads(threads) schedule(static,1)
223
+    for(uint32_t fidx = 0; fidx < folders.size(); fidx++) {
224
+    	string thread_folder = folders[fidx];
225
+
226
+		ParameterList params;
227
+		setDefault(params);	// set default parameters
228
+
229
+		// add input path and output path
230
+		string sequence_path = "", output = "";
231
+		params.file =  path + "/" + thread_folder + "/" + format;
232
+
233
+		vector<int> seq_compression_params;
234
+		seq_compression_params.push_back(CV_IMWRITE_PNG_COMPRESSION);
235
+		seq_compression_params.push_back(0);
236
+		seq_compression_params.push_back(CV_IMWRITE_JPEG_QUALITY);
237
+		seq_compression_params.push_back(100);
238
+
239
+		params.Jets = 1;
240
+
241
+		int start_format = (params.file.find_last_of('/') + 1);
242
+		int end_format = params.file.length() - start_format;
243
+
244
+		sequence_path = params.file.substr(0,start_format);
245
+		string format = params.file.substr(start_format, end_format);
246
+		if(sequence_path[sequence_path.length() - 1] != '/') sequence_path = sequence_path + "/";
247
+
248
+		params.file = sequence_path;
249
+		params.insert("format", format, true);
250
+
251
+		// set output folder!
252
+		output = sequence_path + "/adaptiveFR/";
253
+
254
+		if(sequence_path.empty() || output.empty())
255
+			continue;
256
+
257
+		int len_format = format.find_last_of('.');
258
+		string format_flow = format.substr(0,len_format);
259
+
260
+		params.sequence_start = start;
261
+		for(uint32_t i = 0; i < params.sequence_start_list.size(); i++)
262
+			params.sequence_start_list[i] = start;
263
+
264
+
265
+		if(output[output.length() - 1] != '/') output = output + "/";
266
+
267
+		// set standard stinel params for epic flow
268
+		epic_params_t epic_params;
269
+		epic_params_default(&epic_params);
270
+		variational_params_t flow_params;
271
+		variational_params_default(&flow_params);
272
+		epic_params.pref_nn= 25;
273
+		epic_params.nn= 160;
274
+		epic_params.coef_kernel = 1.1f;
275
+		flow_params.niter_outer = 5;
276
+		flow_params.alpha = 1.0f;
277
+		flow_params.gamma = 0.72f;
278
+		flow_params.delta = 0.0f;
279
+		flow_params.sigma = 1.1f;
280
+
281
+		// create results folder
282
+		boost::filesystem::create_directories(output);
283
+		boost::filesystem::create_directories(output+"tmp/");                // ParameterList result folder
284
+		boost::filesystem::create_directories(output+"sequence/");                // ParameterList result folder
285
+
286
+
287
+		image_t **wx = new image_t*[samples],
288
+				**wy = new image_t*[samples];
289
+		// TODO: TAKE SAMPLES FROM DIFFERENT STEPS
290
+		for(int it = 0; it < samples; it++) {
291
+			if(it > 0) {
292
+				params.sequence_start += params.Jets * sample_step;
293
+				for(uint32_t i = 0; i < params.sequence_start_list.size(); i++) {
294
+					params.sequence_start_list[i] = params.sequence_start; 	//
295
+				}
296
+			}
297
+
298
+			wx[it] = NULL;
299
+			wy[it] = NULL;
300
+			/*
301
+			 * ################### read in image sequence ###################
302
+			 */
303
+			vector<int> red_loc = params.splitParameter<int>("raw_red_loc","1,0");
304
+
305
+			char** img_files = new char*[all_frames];
306
+			color_image_t **seq = new color_image_t*[all_frames];
307
+
308
+			bool success = true;
309
+
310
+			for(uint32_t f = 0; f < all_frames; f++) {
311
+				char img_file[200];
312
+
313
+				if(!sintel) {
314
+					sprintf(img_file, (sequence_path+format).c_str(), params.sequence_start + f * skip);
315
+				} else {
316
+					int sintel_frame = params.sequence_start / 1000;
317
+					int hfr_frame = f * skip + (params.sequence_start % 1000);
318
+
319
+					while(hfr_frame < 0) {
320
+						sintel_frame--;
321
+						hfr_frame = 42 + hfr_frame;
322
+					}
323
+					while(hfr_frame > 41) {
324
+						sintel_frame++;
325
+						hfr_frame = hfr_frame - 42;
326
+					}
327
+
328
+					sprintf(img_file, (sequence_path+format).c_str(), sintel_frame, hfr_frame);
329
+				}
330
+
331
+				if(access(img_file, F_OK) == -1) {
332
+					cerr << "Could not find " << img_file << "!" << endl;
333
+					success = false;
334
+					break;
335
+				}
336
+
337
+				cout << "Reading " << img_file << "..." << endl;
338
+
339
+				Mat img = imread(string(img_file), CV_LOAD_IMAGE_UNCHANGED);         // load images
340
+
341
+				float norm = 1;
342
+				if(img.type() == 2 || params.parameter<bool>("16bit", "0")) {
343
+					norm = 1.0f/255;		// for 16 bit images
344
+					params.insert("16bit", "1", true);
345
+				}
346
+
347
+				// convert to floating point
348
+				img.convertTo(img, CV_32FC(img.channels()));
349
+
350
+				/*
351
+				 * DEMOSAICING
352
+				 */
353
+				if(raw) {
354
+					Mat tmp = img.clone();
355
+					color_image_t* tmp_in = color_image_new(img.cols, img.rows);
356
+					color_image_t* tmp_out = color_image_new(img.cols, img.rows);
357
+
358
+					switch(params.parameter<int>("raw_demosaicing", "0")) {
359
+						case 0: // use bilinear demosaicing
360
+								img = Mat::zeros(tmp.rows, tmp.cols, CV_32FC3);
361
+								bayer2rgbGR(tmp, img, red_loc[0], red_loc[1]); // red green
362
+								break;
363
+
364
+						case 1:	// use hamilton adams demosaicing
365
+								mat2colorImg<float>(img, tmp_in);
366
+
367
+								HADemosaicing(tmp_out->c1, tmp_in->c1, tmp_in->width, tmp_in->height, red_loc[0], red_loc[1]); // Hamilton-Adams implemented by Pascal Getreuer
368
+
369
+								img = Mat::zeros(img.rows, img.cols, CV_32FC3);
370
+								colorImg2colorMat<Vec3f>(tmp_out, img);
371
+								break;
372
+
373
+						case 2: // use opencv demosaicing
374
+								tmp.convertTo(tmp, CV_8UC1);
375
+								img = Mat::zeros(tmp.rows, tmp.cols, CV_8UC3);
376
+
377
+								int code = CV_BayerBG2RGB;
378
+								if(red_loc[1] == 0) // y
379
+									if(red_loc[0] == 0) // x
380
+										code = CV_BayerBG2RGB;
381
+									else
382
+										code = CV_BayerGB2RGB;
383
+								else
384
+									if(red_loc[0] == 0) // x
385
+										code = CV_BayerGR2RGB;
386
+									else
387
+										code = CV_BayerRG2RGB;
388
+
389
+								cv::cvtColor(tmp, img, code); // components from second row, second column !!!!!!!!!!!!!!!!!
390
+								img.convertTo(img, CV_32FC(img.channels()));
391
+								break;
392
+					}
393
+
394
+					color_image_delete(tmp_in);
395
+					color_image_delete(tmp_out);
396
+				} else {
397
+					// covert to RGB
398
+					cv::cvtColor(img, img, CV_BGR2RGB);
399
+				}
400
+
401
+				// use only a part of the images
402
+				if(params.extent.x > 0 || params.extent.y > 0) {
403
+					img = img.rowRange(Range(params.center.y - params.extent.y/2,params.center.y + params.extent.y/2));
404
+					img = img.colRange(Range(params.center.x - params.extent.x/2,params.center.x + params.extent.x/2));
405
+				}
406
+
407
+				// rescale image with gaussian blur to avoid anti-aliasing
408
+				if(scale != 1) {
409
+					GaussianBlur(img, img,Size(),1/sqrt(2*scale),1/sqrt(2*scale),BORDER_REPLICATE);
410
+					resize(img, img, Size(0,0), scale, scale, INTER_LINEAR);
411
+				}
412
+
413
+				// print to file
414
+				img_files[f] = new char[500];
415
+				sprintf(img_files[f], (output+"sequence/frame_%i.png").c_str(), params.sequence_start + f * skip);
416
+				Mat output_img;
417
+
418
+
419
+				if(params.verbosity(WRITE_FILES)) {
420
+					if(params.parameter<bool>("16bit", "0")) {
421
+						img.convertTo(output_img, CV_16UC(img.channels()));
422
+					} else {
423
+						img.convertTo(output_img, CV_8UC(img.channels()), norm);
424
+					}
425
+					cv::cvtColor(output_img, output_img, CV_RGB2BGR);	// OpenCV uses BGR
426
+					imwrite(img_files[f], output_img, seq_compression_params);
427
+				}
428
+
429
+				// use 8 bit for further processing
430
+				img.convertTo(img, CV_8UC(img.channels()), norm);
431
+
432
+				int width = img.cols;
433
+				int height = img.rows;
434
+
435
+				// copy data
436
+				seq[f] = color_image_new(width, height);
437
+
438
+				if(img.channels() == 1) {
439
+					mat2colorImg<uchar>(img, seq[f]);
440
+				} else
441
+					colorMat2colorImg<Vec3b>(img, seq[f]);
442
+
443
+				// resize and copy data for deep match
444
+				GaussianBlur(img, img,Size(),1/sqrt(2*dm_scale),1/sqrt(2*dm_scale),BORDER_REPLICATE);
445
+				resize(img, img, Size(0,0), dm_scale, dm_scale, INTER_LINEAR);
446
+
447
+				sprintf(img_files[f], (output+"sequence/frame_epic_%i.png").c_str(), params.sequence_start + f * skip);
448
+				output_img = Mat(img.rows, img.cols, CV_8UC(img.channels()));
449
+				cv::cvtColor(img, output_img, CV_RGB2BGR);	// OpenCV uses BGR
450
+				imwrite(img_files[f], output_img, seq_compression_params);
451
+			}
452
+
453
+			if(!success)
454
+				continue;
455
+
456
+			/*
457
+			 *  write infos to file
458
+			 */
459
+			params.print();
460
+
461
+			ofstream infos;
462
+			infos.open((output + "config.cfg").c_str());
463
+			infos << "# Epic Flow estimation\n";
464
+			infos << params;
465
+			infos.close();
466
+
467
+			// write stats
468
+			stringstream results;
469
+			results << "frame\ttime\n\n";
470
+			int avg_time = 0;
471
+			int counter = 0;
472
+
473
+			for(uint32_t j = 0; j < params.Jets; j++) {
474
+				ParameterList thread_params(params);
475
+
476
+				int f = j;
477
+
478
+				color_image_t **im = &seq[f];
479
+
480
+				// prepare variables
481
+				wx[it] = image_new(im[0]->width*dm_scale, im[0]->height*dm_scale);
482
+				wy[it] = image_new(im[0]->width*dm_scale, im[0]->height*dm_scale);
483
+
484
+				time_t pp_start, pp_stop;
485
+				int t_preprocessing = 0;
486
+				char edges_f[1000], edges_cmd[1000], match_f[1000], match_cmd[1000], epic_cmd[1000], epic_f[1000];
487
+
488
+				char forward_flow_file[200];
489
+				if(!sintel)
490
+					sprintf(forward_flow_file, (output + format_flow + ".flo").c_str(), params.sequence_start + f * skip);
491
+				else
492
+					sprintf(forward_flow_file, (output + format_flow + ".flo").c_str(), params.sequence_start + f * skip, 0);
493
+
494
+				time_t t_start, t_stop;
495
+				// skip finished frames
496
+				if(overwrite || access( forward_flow_file, F_OK ) == -1) {
497
+					/*
498
+					 * ################### extract edges and get matches ###################
499
+					 */
500
+					cout << "Computing edges ..." << endl;
501
+					sprintf(edges_f, (output+"tmp/edges_%i.dat").c_str(), params.sequence_start + f);
502
+
503
+					if(overwrite || access( edges_f, F_OK ) == -1) {
504
+						sprintf(edges_cmd, "matlab -nodesktop -nojvm -r \"addpath(\'%s/matlab/\'); detect_edges(\'%s\',\'%s\'); exit\"", SOURCE_PATH.c_str(), img_files[j], edges_f);
505
+
506
+						time(&pp_start);
507
+						system(edges_cmd);
508
+						// Call the function
509
+						time(&pp_stop);
510
+						t_preprocessing += (int) difftime(pp_stop, pp_start);
511
+					}
512
+
513
+					// matches to target frame
514
+					cout << "Computing matches between " << params.sequence_start + f * skip << " and " << params.sequence_start + (f + 1) * skip<< " ..." << endl;
515
+					sprintf(match_f, (output+"tmp/matches_%i_%i.dat").c_str(), params.sequence_start + f * skip, params.sequence_start + (f + 1) * skip);
516
+
517
+					cout << img_files[j] << " and " << img_files[j + 1] << endl;
518
+					if(overwrite || access( match_f, F_OK ) == -1) {
519
+						sprintf(match_cmd, "%s/deepmatching %s %s -png_settings -out %s", DEEPMATCHING_PATH.c_str(), img_files[j], img_files[j + 1], match_f);
520
+
521
+						time(&pp_start);
522
+						system(match_cmd);
523
+						time(&pp_stop);
524
+						t_preprocessing += (int) difftime(pp_stop, pp_start);
525
+					}
526
+
527
+
528
+					/*
529
+					 * ############ forward flow ##################
530
+					 */
531
+					cout << "Forward flow estimation ..." << endl;
532
+
533
+					// matches to target frame
534
+					float_image forward_edges = read_edges(edges_f, im[0]->width, im[0]->height);
535
+					time(&pp_start);
536
+					float_image forward_matches = read_matches(match_f);
537
+					time(&pp_stop);
538
+					t_preprocessing += (int) difftime(pp_stop, pp_start);
539
+
540
+					color_image_t *imlab = rgb_to_lab(im[0]);
541
+
542
+					// initilize with deep matches
543
+					cout << "Epic interpolation of forward flow ..." << endl;
544
+					time(&pp_start);
545
+					epic(wx[it], wy[it], imlab, &forward_matches, &forward_edges, &epic_params, 1);
546
+					time(&pp_stop);
547
+					t_preprocessing += (int) difftime(pp_stop, pp_start);
548
+
549
+					// energy minimization
550
+					time(&t_start);
551
+					variational(wx[it], wy[it], im[0], im[1], &flow_params);
552
+					system(epic_cmd);
553
+					time(&t_stop);
554
+					t_preprocessing += difftime(t_stop, t_start);
555
+
556
+					color_image_delete(imlab);
557
+
558
+					free(forward_matches.pixels);
559
+					free(forward_edges.pixels);
560
+
561
+				    // write output file
562
+				    writeFlowFile(forward_flow_file, wx[it], wy[it]);
563
+
564
+					cout << "Forward flow from frame " << params.sequence_start + f * skip << " to " << params.sequence_start + (f + 1) * skip << " finished! (Computation took " << t_preprocessing<< " s)" << endl;
565
+				} else {
566
+					image_t **tmp;
567
+					tmp = readFlowFile(forward_flow_file);
568
+
569
+					wx[it] = tmp[0];
570
+					wy[it] = tmp[1];
571
+
572
+					cout << "Forward flow from frame " << params.sequence_start + f * skip << " to " << params.sequence_start + (f + 1) * skip << " already exist!" << endl;
573
+				}
574
+
575
+				Mat floImg = flowColorImg(wx[it], wy[it], params.verbosity(VER_CMD));
576
+				// write flow image to file
577
+				if(!floImg.data) {                 // Check for invalid input
578
+					cout <<  "No forward flow for frame " << params.sequence_start + f * skip << std::endl ;
579
+				} else {
580
+					if(!output.empty()) {
581
+						stringstream flowF;
582
+						flowF << output <<  "tmp/frame_" << params.sequence_start + f * skip << ".png";
583
+
584
+						imwrite((flowF.str()), floImg);
585
+					}
586
+				}
587
+
588
+				// normalize flow to recorded resolution and frame rate
589
+				image_mul_scalar(wx[it], 1.0f / (scale * skip));		// scale flow
590
+				image_mul_scalar(wy[it], 1.0f / (scale * skip));		// scale flow
591
+			}
592
+
593
+
594
+			if(counter > 0) avg_time /= counter;
595
+
596
+			cout << "Average computation was " << avg_time << " s" << endl;
597
+			results << "\n\navg\t" << avg_time << "s\n";
598
+
599
+			// write experiment results to file
600
+			if(!params.output.empty() && counter > 0) {
601
+				ofstream infos;
602
+				infos.open((output + "results.info").c_str());
603
+				infos << "Epic Flow Multi Frame\n";
604
+				infos << "\n";
605
+				infos << results.str();
606
+				infos.close();
607
+			}
608
+
609
+			// clean up
610
+			for(uint32_t f = 0; f < all_frames; f++) {
611
+				color_image_delete(seq[f]);
612
+				delete[] img_files[f];
613
+			}
614
+			delete[] img_files;
615
+			delete[] seq;
616
+		}
617
+
618
+		/*
619
+		 * ########################################### compute quantil ##############################################
620
+		 */
621
+		int used = 0;
622
+		vector<double> magnitudes;
623
+		magnitudes.reserve(samples * wx[0]->height * wx[0]->width);
624
+		for(int it = 0; it < samples; it++) {
625
+			if(wx[it] == NULL || wy[it] == NULL) continue;
626
+
627
+			for (int y = 0; y < wx[it]->height; y++) {
628
+				for (int x = 0; x < wx[it]->width; x++) {
629
+					magnitudes.push_back(sqrt(wx[it]->data[y * wx[it]->stride + x] * wx[it]->data[y * wx[it]->stride + x] + wy[it]->data[y * wy[it]->stride + x] * wy[it]->data[y * wy[it]->stride + x]));
630
+				}
631
+			}
632
+
633
+			used++;
634
+		}
635
+		sort(magnitudes.begin(), magnitudes.end());
636
+
637
+		float np = q * magnitudes.size() - 1;
638
+		double quantil = 0;
639
+		if((np < magnitudes.size() - 1) && fmod(np,2.0f) == 0) {
640
+			quantil = 0.5f * (magnitudes[(int) np] + magnitudes[(int) np + 1]);
641
+		} else {
642
+			quantil = (magnitudes[(int) ceil(np)]);
643
+		}
644
+
645
+		double maxq = magnitudes.back();
646
+
647
+		cout << "Quantil: " << quantil << endl;
648
+
649
+		// write experiment results to file
650
+		ofstream infos;
651
+		infos.open((output + "results.info").c_str());
652
+		infos << "Adaptive Frame rate\n";
653
+		infos << "\n";
654
+		infos << "samples	" << used << "\n";
655
+		infos << "sample_step	" << sample_step << "\n";
656
+		infos << "skip	" << skip << "\n";
657
+		infos << q << " quantil	" << quantil << "\n";
658
+		infos << "max	" << maxq << "\n";
659
+		infos.close();
660
+
661
+		#pragma omp critical (overview)
662
+		{
663
+			overview << thread_folder << "\t" << q << " quantil\t" << quantil << "\n";
664
+		}
665
+
666
+		string fname = sequence_path + "quantil.dat";
667
+		if(!append.empty())
668
+			infos.open(append.c_str(), std::ofstream::out | std::ofstream::app);
669
+		else
670
+			infos.open(fname.c_str());
671
+		infos << quantil << "\n";
672
+		infos << maxq << "\n";
673
+		infos.close();
674
+
675
+		for(int it = 0; it < samples; it++) {
676
+			image_delete(wx[it]);
677
+			image_delete(wy[it]);
678
+		}
679
+		delete[] wx;
680
+		delete[] wy;
681
+    }
682
+
683
+	ofstream infos;
684
+	infos.open((path + "results.info").c_str());
685
+	infos << "Adaptive Frame rate\n";
686
+	infos << "\n";
687
+	infos << "samples	" << samples << "\n";
688
+	infos << "sample_step	" << sample_step << "\n";
689
+	infos << "skip	" << skip << "\n\n";
690
+	infos << overview.str();
691
+	infos.close();
692
+
693
+	cout << "Done!" << endl;
694
+    return 0;
695
+}

+ 3 - 0
adaptiveFR.dat

@@ -0,0 +1,3 @@
1
+opt_hfr_quantil	2
2
+opt_lfr_quantil	8
3
+opt_lfr_rate	4

+ 58 - 0
cfgs/dense_tracking.cfg

@@ -0,0 +1,58 @@
1
+# Parameters for dense tracking formulation
2
+verbose		00000		# 1. digit - console
3
+threads		1		# specify number of threads (OPENMP necessary)
4
+
5
+file		[path to img sequence]/0%06i.tif	# path and format of image sequence
6
+jet_estimation	[path to slow flow output with adaptive framerate]/low_fr/	# using adaptive frame rate specify path to the low frame rate slow flow estimates
7
+jet_estimation	[path to slow flow output with adaptive framerate]/high_fr/	# using adaptive frame rate specify path to the high frame rate slow flow estimates
8
+flow_format	%07i
9
+
10
+output		[output path]				# output path
11
+
12
+grayscale			0	# set to one if input are grayscale images
13
+16bit				1	# set to 1 if input images 16 bit
14
+raw				1	# set to 1 if raw images			
15
+raw_demosaicing			2	# 0: bilinear green ratio, 1: hamilton adams, 2: opencv
16
+raw_red_loc			1,0	# location of first red pixel (x,y)
17
+
18
+start				10	# first frame in input sequence
19
+adaptive			1	# choose two frame rates according to 0.99 quantile motion
20
+ref_fps_F			450	# number of frames estimated in the final frame rate
21
+max_fps				200	# frame rate of high frame rate sequence
22
+
23
+scale				1.0	# scale input images
24
+acc_skip_pixel			1	# skip number of pixel for final resolution
25
+acc_occlusion			0	# set to 1 to use occlusions from jet optimization
26
+
27
+acc_discard_inconsistent	0	# set to 1 to discard inconsistent trajectories
28
+acc_epic_interpolation		1	# set to 1 to use epic flow interpolation
29
+acc_epic_skip			2	# skip pixel for epic flow interpolation
30
+
31
+acc_jet_consistency		1.0	# weight of flow data term
32
+acc_brightness_constancy	0.1	# weight of brightness constancy
33
+acc_gradient_constancy		1.0	# weight of gradient constancy
34
+acc_occlusion_penalty		500.0	# penalty for occluded frames
35
+acc_beta			10.0	# spatial smoothness flow
36
+acc_spatial_occ			10.0	# spatial smoothness occlusion
37
+acc_cv				0.0	# temporal flow smoothness
38
+acc_temporal_occ		10.0	# temporal occlusion smoothness
39
+
40
+acc_occlusion_threshold		5.0	# occlusion estimation: threshold for comparison of high speed flow and trajectory
41
+acc_occlusion_fb_threshold	5.0	# occlusion estimation: threshold for forward backward consistency check
42
+
43
+acc_alternate			5	# alternate between trws and propagation
44
+acc_approach			0	# optimization 0: TRW-S, 1: BP
45
+acc_trws_eps			1e-5	# stopping criterion for lower bound
46
+acc_trws_max_iter		10	# maximum number of iterations
47
+
48
+acc_neigh_hyp			5	# number of proposals propagated to neighbors
49
+acc_neigh_hyp_radius		100.0	# radius for proposal propagation
50
+acc_hyp_neigh_tryouts		20	# maximum number of tryouts for sampling
51
+
52
+acc_traj_sim_method		1	# compare 0: adjacent flow, 1: accumulated flow, 2: final flow
53
+acc_traj_sim_thres		0.1	# threshold for distance between two hypotheses
54
+
55
+acc_consistency_threshold	1.0	# forward backward consistency threshold
56
+
57
+acc_penalty_fct_data		1	# penality fct data term 0: quadratic 1: modified l1 norm 2: lorentzian
58
+acc_penalty_fct_reg		1	# penality fct regularization term 0: quadratic 1: modified l1 norm 2: lorentzian

+ 60 - 0
cfgs/slow_flow.cfg

@@ -0,0 +1,60 @@
1
+# Default Parameters for Jet Estimation	
2
+verbose		0000				# 1 number - console, 2 - sequence + flow gt, 3 - img pyramid, 4 - flow pyramid
3
+threads		1				# specify number of threads (OPENMP necessary)
4
+
5
+file		[path to img sequence]/0%06i.tif	# path and format of image sequence
6
+
7
+output		[output path]				# output path
8
+
9
+Jets		225				# number of high speed flow estimates
10
+start		10				# first frame of input sequence
11
+max_fps		200				# frame rate of input sequence
12
+ref_fps		20				# final frame rate
13
+
14
+16bit		1				# set to 1 if input images 16 bit
15
+raw		1				# set to 1 if raw images
16
+raw_demosaicing	2				# demosaicing method 0: bilinear interp, 1: hamilton adams, 2: opencv
17
+raw_red_loc	1,0				# location of first red pixel (x,y)
18
+
19
+scale		0.25				# scaling factor for slow_flow
20
+sigma		0				# presmoothing
21
+deep_matching	0				# set to 1 to use deep matching
22
+dm_scale	1.0				# scaling factor for deep matching
23
+
24
+slow_flow_method	symmetric		# symmetric: symmetric window, forward: forward window
25
+slow_flow_S	3				# number of frames in window
26
+
27
+slow_flow_dataterm	1			# 0: unnormalized, 1: normalized
28
+slow_flow_smoothing	1			# 0: \phi(u_dx) + \phi(u_dy), 1: \phi(u_dx + u_symdy) + \phi(u_dy + u_symdx), 2: \phi(u_dx + u_dy)
29
+slow_flow_delta	1.0				# color constancy assumption weight
30
+slow_flow_gamma	6.0				# gradient constancy assumption 
31
+slow_flow_alpha	4.0				# flow smoothness weight
32
+
33
+slow_flow_rho_0	1				# weight for successive data terms
34
+slow_flow_rho_1	1				# weight for successive data terms
35
+slow_flow_omega_0		0		# weight for reference data terms
36
+slow_flow_omega_1		2		# weight for reference data terms
37
+
38
+slow_flow_layers	5			# number of pyramid layers
39
+slow_flow_p_scale	0.9			# scaling factor for pyramid
40
+
41
+slow_flow_occlusion_reasoning	1		# set to 1 to enable occlusion reasoning
42
+slow_flow_occlusion_penalty	0.1		# preference of backwards occlusion (using the forward data terms)
43
+slow_flow_occlusion_alpha	0.1		# occlusion smoothness weight
44
+slow_flow_output_occlusions	1		# set to 1 to output occlusion estimate
45
+
46
+slow_flow_niter_alter	10			# number of alternations 
47
+slow_flow_niter_graphc	10			# number of iterations for graph cut expansion algorithm
48
+slow_flow_niter_outer	10			# number of outer fixed point iterations
49
+slow_flow_thres_outer	1e-5			# threshold for du and dv to stop the optimization
50
+slow_flow_niter_inner	1			# number of inner fixed point iterations
51
+slow_flow_thres_inner	1e-5			# threshold for du and dv to stop the optimization
52
+slow_flow_niter_solver	30			# number of solver iterations
53
+slow_flow_sor_omega	1.9			# omega parameter of sor method
54
+
55
+slow_flow_robust_color	1			# 0: quadratic, 1: modified_l1_norm, 2: lorentzian, 3: trunc mod l1, 4: Geman McClure
56
+slow_flow_robust_color_eps	0.001		# epsilon of robust function
57
+slow_flow_robust_color_truncation	0.5	# truncation threshold (trun_mod_l1)
58
+slow_flow_robust_reg	1			# 0: quadratic, 1: modified_l1_norm, 2: lorentzian, 3: trunc mod l1, 4: Geman McClure
59
+slow_flow_robust_reg_eps	0.001		# epsilon of robust function
60
+slow_flow_robust_reg_truncation	0.5		# truncation threshold (trun_mod_l1)

+ 41 - 0
configuration.h

@@ -0,0 +1,41 @@
1
+/*
2
+ * configuration.h
3
+ *
4
+ *  Created on: Jul 15, 2017
5
+ *      Author: jjanai
6
+ */
7
+
8
+#ifndef CONFIGURATION_H_
9
+#define CONFIGURATION_H_
10
+
11
+#include <iostream>
12
+#include <string>
13
+
14
+// include Hamilton-Adams demosaicing
15
+extern "C"
16
+{
17
+//	#include "[SPECIFY PATH TO HAMILTON ADAMS DEMOSAICING]/dmha.h"
18
+}
19
+
20
+// include flowcode (middlebury devkit)
21
+#include "[SPECIFY PATH TO MIDDLEBURY DEVKIT]/cpp/colorcode.h"
22
+#include "[SPECIFY PATH TO MIDDLEBURY DEVKIT]/cpp/flowIO.h"
23
+
24
+// include TRWS
25
+#include "[SPECIFY PATH TO TRWS]/MRFEnergy.h"
26
+
27
+// path to deepmatching
28
+const string DEEPMATCHING_PATH = "[SPECIFY PATH TO DEEPMATCHING]/deepmatching/";
29
+
30
+// source path
31
+const string SOURCE_FILE = __FILE__;
32
+const string SOURCE_PATH = SOURCE_FILE.substr(0, SOURCE_FILE.rfind("/") + 1);
33
+
34
+void HADemosaicing(float *Output, const float *Input, int Width, int Height, int RedX, int RedY) {
35
+//	HamiltonAdamsDemosaic(Output, Input, Width, Height, RedX, RedY); // Hamilton-Adams implemented by Pascal Getreuer
36
+}
37
+
38
+
39
+
40
+
41
+#endif /* CONFIGURATION_H_ */

+ 15 - 0
configuration_epic.h

@@ -0,0 +1,15 @@
1
+/*
2
+ * configuration_epic.h
3
+ *
4
+ *  Created on: Jul 15, 2017
5
+ *      Author: jjanai
6
+ */
7
+
8
+#ifndef CONFIGURATION_EPIC_H_
9
+#define CONFIGURATION_EPIC_H_
10
+
11
+// include graph cut library
12
+#include "[SPECIFY PATH TO GCO]/GCoptimization.h"
13
+
14
+
15
+#endif /* CONFIGURATION_EPIC_H_ */

+ 1961 - 0
dense_tracking.cpp

@@ -0,0 +1,1961 @@
1
+/*
2
+ * dense_tracking.cpp
3
+ *
4
+ *  Created on: Mar 7, 2016
5
+ *      Author: Janai
6
+ */
7
+#include <iostream>
8
+#include <fstream>
9
+#include <stdio.h>
10
+#include <string>
11
+#include <cmath>
12
+#include <omp.h>
13
+#include <random>
14
+#include <set>
15
+
16
+#include <flann/flann.hpp>
17
+
18
+extern "C" {
19
+	#include "../libs/dmgunturk/dmha.h"
20
+}
21
+#include <chrono>
22
+
23
+#include <opencv2/core.hpp>
24
+#include <opencv2/highgui.hpp>
25
+#include <opencv2/imgproc.hpp>
26
+
27
+#include <boost/filesystem.hpp>
28
+#include <boost/multi_index_container.hpp>
29
+#include <boost/multi_index/indexed_by.hpp>
30
+#include <boost/multi_index/ordered_index.hpp>
31
+#include <boost/multi_index/composite_key.hpp>
32
+#include <boost/multi_index/member.hpp>
33
+#include <boost/chrono/thread_clock.hpp>
34
+
35
+#include <gsl/gsl_multifit.h>
36
+#include <gsl/gsl_randist.h>
37
+
38
+#include "epic_flow_extended/image.h"
39
+#include "epic_flow_extended/io.h"
40
+#include "epic_flow_extended/epic.h"
41
+#include "epic_flow_extended/variational_mt.h"
42
+#include "utils/utils.h"
43
+#include "utils/hypothesis.h"
44
+#include "utils/parameter_list.h"
45
+#include "penalty_functions/penalty_function.h"
46
+#include "penalty_functions/lorentzian.h"
47
+#include "penalty_functions/modified_l1_norm.h"
48
+#include "penalty_functions/quadratic_function.h"
49
+#include "configuration.h"
50
+
51
+
52
+using namespace std;
53
+using namespace cv;
54
+using namespace boost::multi_index;
55
+using namespace boost::chrono;
56
+
57
+void setDefaultVariational(ParameterList& params) {
58
+    // general
59
+    params.insert("verbose", "0", true);
60
+    params.insert("threads", "1", true);
61
+
62
+    params.insert("scale", "1.0f", true);
63
+
64
+    params.insert("slow_flow_S", "2", true);
65
+
66
+    // energy function
67
+    params.insert("slow_flow_alpha", "4.0f");
68
+    params.insert("slow_flow_gamma", "6.0f", true);
69
+    params.insert("slow_flow_delta", "1.0f", true);
70
+
71
+    // image pyramid
72
+    params.insert("slow_flow_layers", "1", true);
73
+    params.insert("slow_flow_p_scale", "0.9f", true);
74
+
75
+    // optimization
76
+    params.insert("slow_flow_niter_alter", "10", true);
77
+    params.insert("slow_flow_niter_outer", "10", true);
78
+    params.insert("slow_flow_thres_outer", "1e-5", true);
79
+    params.insert("slow_flow_niter_inner", "1", true);
80
+    params.insert("slow_flow_thres_inner", "1e-5", true);
81
+    params.insert("slow_flow_niter_solver", "30", true);
82
+    params.insert("slow_flow_sor_omega", "1.9f", true);
83
+
84
+    // occlusion reasoning
85
+    params.insert("slow_flow_occlusion_reasoning", "1", true);
86
+    params.insert("slow_flow_occlusion_penalty", "0.1", true);
87
+    params.insert("slow_flow_occlusion_alpha", "0.1", true);
88
+    params.insert("slow_flow_output_occlusions", "1", true);
89
+
90
+    // regularization
91
+    params.insert("slow_flow_robust_color", "1", true);
92
+    params.insert("slow_flow_robust_color_eps", "0.001", true);
93
+    params.insert("slow_flow_robust_color_truncation", "0.5", true);
94
+    params.insert("slow_flow_robust_reg", "1", true);
95
+    params.insert("slow_flow_robust_reg_eps", "0.001", true);
96
+    params.insert("slow_flow_robust_reg_truncation", "0.5", true);
97
+}
98
+
99
+void setDefault(ParameterList& params) {
100
+    params.insert("verbose", "0", true);
101
+	// dense tracking parameters
102
+	params.insert("scale", "1");						// scale input images
103
+	params.insert("acc_skip_pixel", "1");				// skip number of pixel for final resolution)
104
+	params.insert("acc_occlusion", "0");				// set to 1 to use occlusions from jet optimization
105
+
106
+	// data driven proposal generation
107
+	params.insert("acc_consistency_threshold", "1.0");  // threshold for forward backward consistency
108
+	params.insert("acc_discard_inconsistent", "1");		// set to 1 to discard inconsistent trajectories
109
+	params.insert("acc_epic_interpolation", "1");		// set to 1 to use epic flow interpolation
110
+	params.insert("acc_epic_skip", "2");				// skip pixel for epic flow interpolation
111
+
112
+	// energy weights and penalites
113
+	params.insert("acc_jet_consistency", "1.0");		// weight of flow data term
114
+	params.insert("acc_brightness_constancy", "0.1");	// weight of brightness constancy
115
+	params.insert("acc_gradient_constancy", "1.0");		// weight of gradient constancy
116
+	params.insert("acc_occlusion_penalty", "500.0");	// penalty for occluded frames
117
+	params.insert("acc_beta", "10.0");					// spatial smoothness flow
118
+	params.insert("acc_satial_occ", "10.0");			// spatial smoothness occlusion
119
+	params.insert("acc_temporal_occ", "10.0");			// temporal smoothness occlusion
120
+	params.insert("acc_cv", "0.0");						// temporal flow smoothness
121
+
122
+	params.insert("acc_traj_sim_method", "1");			// spatial trajectory smoothness compares 0: adjacent flow, 1: accumulated flow, 2: final flow
123
+	params.insert("acc_traj_sim_thres", "0.1");			// threshold for distance between two hypotheses
124
+
125
+	// occlusion initialization
126
+	params.insert("acc_occlusion_threshold", "5.0");	// threshold for comparison of high speed flow and trajectory
127
+	params.insert("acc_occlusion_fb_threshold", "5.0");	// threshold for forward backward consistency check
128
+
129
+	// discrete optimization
130
+	params.insert("acc_alternate", "5");				// alternate between trws and propagation
131
+	params.insert("acc_approach", "0");					// 0: TRW-S, 1: BP
132
+	params.insert("acc_trws_eps", "1e-5");				// stopping criterion for lower bound
133
+	params.insert("acc_trws_max_iter", "10");			// maximum number of iterations
134
+
135
+	// neighbor propagation
136
+	params.insert("acc_neigh_hyp", "5");				// number of proposals propagated to neighbors
137
+	params.insert("acc_neigh_hyp_radius", "100.0");		// radius for proposal propagation
138
+	params.insert("acc_neigh_skip1", "2");				// skipped pixel for nearest neighbor search tree
139
+	params.insert("acc_neigh_skip2", "4");				// skipped pixel for nearest neighbor search tree
140
+	params.insert("acc_hyp_neigh_tryouts", "20");		// maximum number of tryouts for sampling
141
+
142
+	// penalty
143
+	params.insert("acc_penalty_fct_data", "1");			// penality fct data term 0: quadratic 1: modified l1 norm 2: lorentzian
144
+	params.insert("acc_penalty_fct_data_eps", "0.001");
145
+	params.insert("acc_penalty_fct_reg", "1");			// penality fct regularization term 0: quadratic 1: modified l1 norm 2: lorentzian
146
+	params.insert("acc_penalty_fct_reg_eps", "0.001");
147
+}
148
+
149
+inline bool insideImg(double x, double y, int width, int height) {
150
+	return (y >= 0 && y < height && x >= 0 && x < width);
151
+}
152
+
153
+bool compareHypotheses(hypothesis* h1, hypothesis* h2) {
154
+	return (h1 != NULL && (h2 == NULL || h1->score() < h2->score()));	// ascending order considering the score and null hypothesis to the end
155
+}
156
+
157
+float addJC(hypothesis* h, const Mat* obs, double acc_jc, double acc_cv, PenaltyFunction* phi_d, ParameterList& params, const Mat* occlusion_masks = NULL) {
158
+	Point2d p = h->p;
159
+
160
+	int height = obs[0].rows;
161
+	int width = obs[0].cols;
162
+
163
+	double jenergy = 0;
164
+	double cvenergy = 0;
165
+	int contribution = 0;
166
+
167
+	for (uint32_t j = 0; j < params.Jets; j++) {
168
+		double u_j = h->u(j);
169
+		double v_j = h->v(j);  // x, y
170
+
171
+		double u_jm1 = 0;
172
+		double v_jm1 = 0;
173
+		if (j > 0) {
174
+			u_jm1 = h->u((j - 1));
175
+			v_jm1 = h->v((j - 1));  // x, y
176
+		}
177
+
178
+		// exclude unknown flow
179
+		if (u_j > UNKNOWN_FLOW_THRESH || v_j > UNKNOWN_FLOW_THRESH)
180
+			break;
181
+
182
+		// ############################################### compare flow to jet estimation
183
+		if (insideImg(p.x + u_jm1, p.y + v_jm1, width, height)) {
184
+			if (h->occluded(j) == 1 || h->occluded(j+1) == 1)	// flow from j to j+1
185
+				continue;
186
+
187
+			double I_x = bilinearInterp<double>(p.x + u_jm1, p.y + v_jm1, obs[j], 1);
188
+			double I_y = bilinearInterp<double>(p.x + u_jm1, p.y + v_jm1, obs[j], 0);  // x, y
189
+
190
+			jenergy += 0.5 * phi_d->apply(((u_j - u_jm1 - I_x) * (u_j - u_jm1 - I_x) + (v_j - v_jm1 - I_y) * (v_j - v_jm1 - I_y)));
191
+
192
+			contribution++;
193
+		}
194
+
195
+		double u_jp1 = 0;
196
+		double v_jp1 = 0;
197
+		if (j + 1 < params.Jets) {
198
+			u_jp1 = h->u((j + 1));
199
+			v_jp1 = h->v((j + 1));  // x, y
200
+		}
201
+
202
+		// ############################################### assume constant velocity
203
+		double u_sq = 2 * u_j - u_jm1 - u_jp1;
204
+		double v_sq = 2 * v_j - v_jm1 - v_jp1;
205
+		u_sq *= u_sq;
206
+		v_sq *= v_sq;
207
+
208
+		cvenergy += sqrt(u_sq + v_sq);
209
+	}
210
+
211
+	if(contribution > 0) jenergy /= contribution;
212
+
213
+	return acc_jc * jenergy + acc_cv * cvenergy;
214
+}
215
+
216
+double dt_warp_time = 0, dt_med_time = 0, dt_sum_time = 0;
217
+
218
+/*
219
+ *  ################################### CONSIDER ALL POSSIBLE PAIRWISE TERMS #########################################
220
+ */
221
+float addBCGC(hypothesis* h, color_image_t const* const * obs, color_image_t const* const * dx, color_image_t const* const * dy, double acc_bc, double acc_gc, int skip, ParameterList& params, const Mat* occlusion_masks = NULL) {
222
+	Point2d p = h->p;
223
+
224
+	int r = 0.5f * (skip + 1);
225
+
226
+	int height = obs[0]->height;
227
+	int width = obs[0]->width;
228
+	int stride = obs[0]->stride;
229
+
230
+	double wenergy = 0;
231
+	double neighs = 0;
232
+
233
+	time_t dt_start, dt_end;
234
+
235
+	// ############################# Mean brightness and gradient constancy #############################
236
+	for (int off_x = (p.x - r); off_x <= (p.x + r); off_x++) {
237
+		for (int off_y = (p.y - r); off_y <= (p.y + r); off_y++) {
238
+			if (off_x < 0 || off_x >= width || off_y < 0 || off_y >= height)
239
+				continue;
240
+
241
+			// warp images
242
+			uint32_t visible = 0;
243
+			vector<vector<double> > I(3, vector<double>(params.Jets + 1, 0));
244
+			vector<vector<double> > Ix(3, vector<double>(params.Jets + 1, 0));
245
+			vector<vector<double> > Iy(3, vector<double>(params.Jets + 1, 0));
246
+
247
+			time(&dt_start);
248
+			for (uint32_t j = 0; j < params.Jets + 1; j++) {
249
+				double x_j = off_x;
250
+				double y_j = off_y;
251
+
252
+				if (j == 0) {
253
+					int idx = stride * y_j + x_j;
254
+					I[0][j] = obs[j]->c3[idx];
255
+					I[1][j] = obs[j]->c2[idx];
256
+					I[2][j] = obs[j]->c1[idx];
257
+					Ix[0][j] = dx[j]->c3[idx];
258
+					Ix[1][j] = dx[j]->c2[idx];
259
+					Ix[2][j] = dx[j]->c1[idx];
260
+					Iy[0][j] = dy[j]->c3[idx];
261
+					Iy[1][j] = dy[j]->c2[idx];
262
+					Iy[2][j] = dy[j]->c1[idx];
263
+
264
+					visible++;
265
+				} else {
266
+					x_j += h->u(j - 1);
267
+					y_j += h->v(j - 1);
268
+
269
+					// stop in the case the object will stay occluded!
270
+					if (insideImg(x_j, y_j, width, height) && (occlusion_masks == NULL || occlusion_masks[j].at<uchar>(y_j, x_j) != 0)) {
271
+						I[0][j] = bilinearInterp(x_j, y_j, obs[j]->c3, height, width, stride);
272
+						I[1][j] = bilinearInterp(x_j, y_j, obs[j]->c2, height, width, stride);
273
+						I[2][j] = bilinearInterp(x_j, y_j, obs[j]->c1, height, width, stride);
274
+						Ix[0][j] = bilinearInterp(x_j, y_j, dx[j]->c3, height, width, stride);
275
+						Ix[1][j] = bilinearInterp(x_j, y_j, dx[j]->c2, height, width, stride);
276
+						Ix[2][j] = bilinearInterp(x_j, y_j, dx[j]->c1, height, width, stride);
277
+						Iy[0][j] = bilinearInterp(x_j, y_j, dy[j]->c3, height, width, stride);
278
+						Iy[1][j] = bilinearInterp(x_j, y_j, dy[j]->c2, height, width, stride);
279
+						Iy[2][j] = bilinearInterp(x_j, y_j, dy[j]->c1, height, width, stride);
280
+
281
+						visible++;
282
+					}
283
+				}
284
+			}
285
+			time(&dt_end);
286
+			dt_warp_time += difftime(dt_end, dt_start);
287
+
288
+			int contribution = 0;
289
+			double e_p = 0;
290
+
291
+
292
+			time(&dt_start);
293
+			// compute data terms
294
+			for (uint32_t i = 0; i < visible; i++) {
295
+				for (uint32_t j = (i + 1); j < visible; j++) {
296
+					double x_i = off_x;
297
+					double y_i = off_y;
298
+					if (i > 0) {
299
+						x_i += h->u(i - 1);
300
+						y_i += h->v(i - 1);
301
+					}
302
+					double x_j = off_x + h->u(j - 1);
303
+					double y_j = off_y + h->v(j - 1);  // x, y
304
+
305
+					if (insideImg(x_i, y_i, width, height) && insideImg(x_j, y_j, width, height)) {
306
+						if (h->occluded(i) == 1 || h->occluded(j) == 1)
307
+							continue;															// skip occluded jets
308
+
309
+						e_p += acc_bc * 0.3334  * (fabs(I[0][i] - I[0][j]) + fabs(I[1][i] - I[1][j]) + fabs(I[2][i] - I[2][j]));
310
+						e_p += acc_gc * 0.3334 * (fabs(Ix[0][i] - Ix[0][j]) + fabs(Ix[1][i] - Ix[1][j]) + fabs(Ix[2][i] - Ix[2][j]) + fabs(Iy[0][i] - Iy[0][j]) + fabs(Iy[1][i] - Iy[1][j]) + fabs(Iy[2][i] - Iy[2][j]));
311
+
312
+						contribution++;
313
+					}
314
+				}
315
+			}
316
+			time(&dt_end);
317
+			dt_sum_time += difftime(dt_end, dt_start);
318
+
319
+
320
+			if(contribution > 0) e_p /= contribution;
321
+
322
+			wenergy += e_p;
323
+			neighs++;
324
+		}
325
+	}
326
+
327
+	if(neighs > 0) wenergy /= neighs;
328
+
329
+	return wenergy;
330
+}
331
+
332
+float addOC(hypothesis* h, double acc_occ, double acc_temporal_occ, ParameterList& params) {
333
+	int occlusions = 0;
334
+	int change = 0;
335
+
336
+	for (uint32_t i = 0; i < params.Jets + 1; i++) {
337
+		// ############################# prefer smaller number of occlusions #############################
338
+		occlusions += h->occluded(i);															// skip occluded jets
339
+
340
+		// ############################# expect temporal consistancy #############################
341
+		if (i < params.Jets && h->occluded(i) != h->occluded(i + 1))
342
+			change++;
343
+	}
344
+
345
+	return acc_occ * occlusions + acc_temporal_occ * change;
346
+}
347
+
348
+void computeSmoothnessWeight(const color_image_t *im, image_t* lum, float coef, const convolution_t *deriv, float avg_1, float avg_2, float avg_3, float std_1, float std_2, float std_3, bool hbit) {
349
+	int i;
350
+	image_t *lum_x = image_new(im->width, im->height), *lum_y = image_new(im->width, im->height);
351
+
352
+	// compute luminance
353
+	v4sf *im1p = (v4sf*) im->c1, *im2p = (v4sf*) im->c2, *im3p = (v4sf*) im->c3, *lump = (v4sf*) lum->data;
354
+	for (i = 0; i < im->height * im->stride / 4; i++) {
355
+		if (hbit)
356
+			*lump = (0.299f * ((*im1p) * std_1 + avg_1) + 0.587f * ((*im2p) * std_2 + avg_2) + 0.114f * ((*im3p) * std_3 + avg_3)) / 65535.0f;	// channels normalized to 0..1
357
+		else
358
+			*lump = (0.299f * ((*im1p) * std_1 + avg_1) + 0.587f * ((*im2p) * std_2 + avg_2) + 0.114f * ((*im3p) * std_3 + avg_3)) / 255.0f;	// channels normalized to 0..1
359
+
360
+		lump += 1;
361
+		im1p += 1;
362
+		im2p += 1;
363
+		im3p += 1;
364
+	}
365
+
366
+	// compute derivatives with five-point tencil
367
+	convolve_horiz(lum_x, lum, deriv);
368
+	convolve_vert(lum_y, lum, deriv);
369
+
370
+	// compute lum norm
371
+	lump = (v4sf*) lum->data;
372
+	v4sf *lumxp = (v4sf*) lum_x->data, *lumyp = (v4sf*) lum_y->data;
373
+	for (i = 0; i < lum->height * lum->stride / 4; i++) {
374
+		*lump = -coef * __builtin_ia32_sqrtps((*lumxp) * (*lumxp) + (*lumyp) * (*lumyp));
375
+		lump[0][0] = 0.5f * expf((float) lump[0][0]);
376
+		lump[0][1] = 0.5f * expf((float) lump[0][1]);
377
+		lump[0][2] = 0.5f * expf((float) lump[0][2]);
378
+		lump[0][3] = 0.5f * expf((float) lump[0][3]);
379
+
380
+		lump += 1;
381
+		lumxp += 1;
382
+		lumyp += 1;
383
+	}
384
+
385
+	image_delete(lum_x);
386
+	image_delete(lum_y);
387
+}
388
+
389
+/* show usage information */
390
+void usage(){
391
+    printf("usage:\n");
392
+    printf("    ./dense_tracking [cfg] -select [estimation for one specific final pair] -resume\n");
393
+    printf("\n");
394
+}
395
+
396
+int main(int argc, char **argv) {
397
+	if (argc < 2) {
398
+		usage();
399
+		return -1;
400
+	}
401
+
402
+	/*
403
+	 *  getting config file location
404
+	 */
405
+	uint32_t selected = 0;
406
+	uint32_t selected_end = 0;
407
+
408
+	/*
409
+	 *  read in parameters and print config
410
+	 */
411
+	string file = string(argv[1]);
412
+
413
+	if (boost::filesystem::exists(file)) {
414
+		printf("using parameters %s\n", file.c_str());
415
+	} else {
416
+		usage();
417
+		return -1;
418
+	}
419
+
420
+	ParameterList params;
421
+	setDefault(params);
422
+	params.read(file);
423
+
424
+    bool resume_frame = false;
425
+	int max_fps = params.parameter<int>("max_fps", "0");		// fps of original sequence
426
+	int threads = params.parameter<int>("threads", "1");
427
+
428
+	#define isarg(key)  !strcmp(a,key)
429
+	if (argc > 1) {
430
+		int current_arg = 1;
431
+		while (current_arg < argc) {
432
+			const char* a = argv[current_arg++];
433
+			if (a[0] != '-') {
434
+				continue;
435
+			}
436
+
437
+			if ( isarg("-h") || isarg("-help"))
438
+				usage();
439
+			else if (isarg("-output"))
440
+				params.output = string(argv[current_arg++]);
441
+			else if (isarg("-threads"))
442
+				threads = atoi(argv[current_arg++]);
443
+			else if( isarg("-resume") )
444
+				resume_frame = true;
445
+			else if( isarg("-select") ) {
446
+				selected = atoi(argv[current_arg++]);
447
+				selected_end = selected + 1;
448
+			} else {
449
+				fprintf(stderr, "unknown argument %s", a);
450
+				usage();
451
+				exit(1);
452
+			}
453
+		}
454
+	} else {
455
+		usage();
456
+		return -1;
457
+	}
458
+
459
+	// check path
460
+	for (uint32_t i = 0; i < params.jet_estimation.size(); i++)
461
+		if (params.jet_estimation[i].back() != '/') params.jet_estimation[i] += "/";
462
+
463
+
464
+	bool sintel = params.parameter<bool>("sintel", "0");
465
+	bool subframes = params.parameter<bool>("subframes", "0");	// are subframes specified
466
+
467
+	int skip_pixel = params.parameter<int>("acc_skip_pixel", "0");
468
+
469
+	int ref_fps_F = params.parameter<int>("ref_fps_F", "1");
470
+	uint32_t rates = params.jet_estimation.size();
471
+	vector<float> weight_jet_estimation(rates, 0);
472
+	for (uint32_t i = 0; i < rates; i++) {
473
+		weight_jet_estimation[i] = i;
474
+		if(params.jet_weight.size() > i)
475
+			weight_jet_estimation[i] = params.jet_weight[i];
476
+	}
477
+	int min_fps_idx = params.parameter<int>("acc_min_fps", "0");
478
+
479
+	/*
480
+	 * ###############################  get step sizes ###############################
481
+	 */
482
+	if (rates != params.jet_S.size()) {
483
+		bool error = false;
484
+		params.jet_S.resize(rates);
485
+		for(uint32_t r = 0; r < rates; r++) {
486
+			if(access((params.jet_estimation[r] + "/config.cfg").c_str(), F_OK) != -1) {
487
+				ParameterList tmp((params.jet_estimation[r] + "/config.cfg"));
488
+
489
+				if(tmp.exists("slow_flow_S"))
490
+					params.jet_S[r] = tmp.parameter<int>("slow_flow_S");
491
+				else {
492
+					cerr << "Error reading S from " << params.jet_estimation[r] << "/config.cfg" << endl;
493
+					error = true;
494
+				}
495
+
496
+			} else {
497
+				cerr << "Error reading " << params.jet_estimation[r] << "/config.cfg" << endl;
498
+				error = true;
499
+			}
500
+		}
501
+
502
+		if(error) {
503
+			cerr << "Frame rate or window size for Jet estimations missing!" << endl;
504
+			exit(-1);
505
+		}
506
+	}
507
+
508
+	// first slow flow estimation is used as reference
509
+	int steps = params.jet_S[min_fps_idx] - 1;
510
+
511
+	/*
512
+	 *  get frame rates
513
+	 */
514
+	if (rates <= 0) {
515
+		cerr << "No Jet estimation specified!" << endl;
516
+		exit(-1);
517
+	}
518
+	if (rates != params.jet_fps.size()) {
519
+		bool error = false;
520
+		params.jet_fps.resize(rates);
521
+
522
+		for(uint32_t r = 0; r < rates; r++) {
523
+			if(access((params.jet_estimation[r] + "/config.cfg").c_str(), F_OK) != -1) {
524
+				ParameterList tmp((params.jet_estimation[r] + "/config.cfg"));
525
+
526
+				if(tmp.exists("jet_fps")) {
527
+					params.jet_fps[r] = tmp.parameter<int>("jet_fps");
528
+				} else {
529
+					cerr << "Error reading jet_fps from " << params.jet_estimation[r] << "/config.cfg" << endl;
530
+					error = true;
531
+				}
532
+			} else {
533
+				cerr << "Error reading " << params.jet_estimation[r] << "/config.cfg" << endl;
534
+				error = true;
535
+			}
536
+		}
537
+
538
+		if(error) {
539
+			cerr << "Frame rate or window size for Jet estimations missing!" << endl;
540
+			exit(-1);
541
+		}
542
+	}
543
+
544
+	// ############################### specify reference framerate by maximum flow ########################################################
545
+	params.Jets = params.jet_fps[min_fps_idx] / (1.0f * params.parameter<int>("ref_fps") * steps);
546
+
547
+	/*
548
+	 * decompose sequence in subsets and adjust number of frames if necessary
549
+	 */
550
+	uint32_t Jets = params.Jets;
551
+
552
+	int skip = (1.0f * max_fps) / params.jet_fps[min_fps_idx];													// number of frames to skip for jet fps
553
+
554
+	uint32_t gtFrames = params.Jets * skip;
555
+
556
+	vector<int> compression_params;
557
+	compression_params.push_back(CV_IMWRITE_PNG_COMPRESSION);
558
+	compression_params.push_back(0);
559
+	compression_params.push_back(CV_IMWRITE_JPEG_QUALITY);
560
+	compression_params.push_back(100);
561
+
562
+	// MAKE SURE FOLDER IS NOT OVERWRITTEN
563
+	if(!resume_frame) {
564
+		string newPath = params.output;
565
+		if(newPath[newPath.length() - 1] == '/') newPath.erase(newPath.length() - 1);
566
+		int num = 1;
567
+		while(boost::filesystem::exists(newPath)) {
568
+			cerr << newPath << " already exists!" << endl;
569
+			stringstream tmp;
570
+			tmp << newPath << "_" << num++;
571
+			newPath = tmp.str();
572
+		}
573
+		params.output = newPath;
574
+	}
575
+	if(params.output[params.output.length() - 1] != '/') params.output = params.output + "/";
576
+
577
+	boost::filesystem::create_directories(params.output);               				// accumulate folder
578
+
579
+	ofstream infos;
580
+	infos.open((params.output + "/config.cfg").c_str());
581
+	infos << "# Slow Flow Accumulation\n";
582
+	infos << params.cfgString(true);
583
+	infos.close();
584
+	cout << (params.output + "config.cfg").c_str() << endl;
585
+
586
+	double traj_sim_thres = params.parameter<double>("acc_traj_sim_thres");
587
+	int traj_sim_method = params.parameter<int>("acc_traj_sim_method");
588
+
589
+	bool use_oracle = params.parameter<bool>("acc_oracle");
590
+	bool use_occlusion = params.parameter<bool>("acc_occlusion");
591
+	float acc_jc = params.parameter<float>("acc_jet_consistency");
592
+	float acc_bc = params.parameter<float>("acc_brightness_constancy");
593
+	float acc_gc = params.parameter<float>("acc_gradient_constancy");
594
+	float acc_occ = params.parameter<float>("acc_occlusion_penalty");
595
+
596
+	double acc_beta = params.parameter<double>("acc_beta");										// weight for flow spatial smoothness
597
+	double acc_spatial_occ = params.parameter<double>("acc_spatial_occ");								// weight for occlusion spatial smoothness
598
+	double acc_temporal_occ = params.parameter<double>("acc_temporal_occ");								// weight for occlusion spatial smoothness
599
+	double acc_cv = params.parameter<double>("acc_cv");											// weight for constant velocity
600
+	double outlier_beta = params.parameter<double>("acc_outlier_beta");							// energy for occlusion penalty
601
+
602
+	float occlusion_threshold = params.parameter<float>("acc_occlusion_threshold");							// hypotheses from neighbors
603
+	float occlusion_fb_threshold = params.parameter<float>("acc_occlusion_fb_threshold");							// hypotheses from neighbors
604
+
605
+	// propagate hypotheses to neighbors
606
+	int alternate = params.parameter<int>("acc_alternate");									// alternate between trws and perturbation
607
+	uint32_t perturb_keep = params.parameter<int>("acc_perturb_keep");							// top x hypotheses will be kept
608
+	bool discard_inconsistent = params.parameter<bool>("acc_discard_inconsistent");
609
+	bool use_jet_occlusions = params.parameter<bool>("acc_use_jet_occlusions");
610
+
611
+	uint32_t hyp_neigh = params.parameter<int>("acc_neigh_hyp");									// hypotheses from neighbors
612
+	float hyp_neigh_radius = params.parameter<int>("acc_neigh_hyp_radius");					// radius to draw neighbors from
613
+	bool draw_nn_from_radius = (hyp_neigh_radius > 0);
614
+	float hyp_neigh_draws = params.parameter<int>("acc_neigh_draws");					// radius to draw neighbors from
615
+	int hyp_neigh_tryouts = params.parameter<int>("acc_hyp_neigh_tryouts");					// maximum number of tryouts for sampling
616
+
617
+	int nn_skip1 = params.parameter<int>("acc_neigh_skip1", "2");
618
+	int nn_skip2 = params.parameter<int>("acc_neigh_skip2", "4");
619
+
620
+	// discrete optimization robust penalty function
621
+	int penalty_fct_data = params.parameter<int>("acc_penalty_fct_data");						// 0:QUADRATIC, 1:MOD_L1, 2:LORETZIAN
622
+	double penalty_fct_data_eps = params.parameter<double>("acc_penalty_fct_data_eps");
623
+	int penalty_fct_reg = params.parameter<int>("acc_penalty_fct_reg");							// 0:QUADRATIC, 1:MOD_L1, 2:LORETZIAN
624
+	double penalty_fct_reg_eps = params.parameter<double>("acc_penalty_fct_reg_eps");
625
+
626
+	unsigned seed = std::chrono::system_clock::now().time_since_epoch().count();
627
+	if (params.exists("seed"))
628
+		seed = params.parameter<int>("seed");
629
+	default_random_engine generator(seed);
630
+
631
+	// EPIC FLOW
632
+	float epic_skip = params.parameter<float>("acc_epic_skip", "2");
633
+	epic_params_t epic_params;
634
+	epic_params_default(&epic_params);
635
+	epic_params.pref_nn = 25;
636
+	epic_params.nn = 160;
637
+	epic_params.coef_kernel = 1.1f;
638
+
639
+	// discrete optimization set options
640
+	MRFEnergy<TypeGeneral>::Options options;
641
+	options.m_method = params.parameter<int>("acc_approach");  							// 0: TRW-S, 1: BP
642
+	options.m_eps = params.parameter<double>("acc_trws_eps");
643
+	options.m_iterMax = params.parameter<int>("acc_trws_max_iter");
644
+	options.m_verbosityLevel = params.verbosity(VER_CMD);
645
+	options.m_printIter = 1;
646
+	options.m_printMinIter = 0;
647
+
648
+	// nearest neighbor search
649
+	flann::SearchParams flann_params;
650
+	flann_params.checks = 128;				//
651
+	flann_params.max_neighbors = -1;		// unlimited
652
+	flann_params.sorted = true;
653
+
654
+	flann::SearchParams flann_params_p;
655
+	flann_params_p.checks = 128;				//
656
+	flann_params_p.max_neighbors = -1;		// unlimited
657
+	flann_params_p.sorted = true;
658
+
659
+	// discrete optimization set robust function
660
+	PenaltyFunction *phi_d = NULL, *phi_s = NULL;
661
+	switch (penalty_fct_data) {
662
+		case 0:
663
+			phi_d = new QuadraticFunction();
664
+			break;
665
+		case 1:
666
+			phi_d = new ModifiedL1Norm(penalty_fct_data_eps);
667
+			break;
668
+		default:
669
+			phi_d = new Lorentzian(penalty_fct_data_eps);
670
+			break;
671
+	}
672
+	switch (penalty_fct_reg) {
673
+		case 0:
674
+			phi_s = new QuadraticFunction();
675
+			break;
676
+		case 1:
677
+			phi_s = new ModifiedL1Norm(penalty_fct_reg_eps);
678
+			break;
679
+		default:
680
+			phi_s = new Lorentzian(penalty_fct_reg_eps);
681
+			break;
682
+	}
683
+
684
+	stringstream acc_folder;
685
+
686
+	acc_folder << params.output;
687
+
688
+	if(use_oracle)
689
+		cout << endl << "Oracle uses same penalty as first jet estimation!!!!" << endl;
690
+
691
+	boost::filesystem::create_directories(acc_folder.str());
692
+	boost::filesystem::create_directories(acc_folder.str() + "gt_occlusions/"); // accumulate folder
693
+	boost::filesystem::create_directories(acc_folder.str() + "occlusions/");  	// accumulate folder
694
+	boost::filesystem::create_directories(acc_folder.str() + "tmp/");           // accumulate folder
695
+
696
+	params.insert("accumulation", acc_folder.str() + "frame_%i.flo", true);
697
+
698
+	params.print();
699
+
700
+	// in case of sintel data we would like to be able to distinguish frame number from 24 fps and 1008 fps
701
+	if (sintel && !subframes)
702
+		params.sequence_start = params.sequence_start * 1000; 	//
703
+
704
+
705
+	if(selected_end == 0)
706
+		selected_end = ref_fps_F;
707
+
708
+	// iterate while changing start
709
+	#pragma omp parallel for num_threads(threads) schedule(static,1)
710
+	for(uint32_t start_jet = selected; start_jet < selected_end; start_jet++) {
711
+
712
+		// alternate between continuous and discrete optimization
713
+		stringstream numVariablesStream, factorsStream;
714
+		double avg_unary_time = 0, avg_pairw_time = 0, avg_optimization_time = 0;
715
+
716
+		// thread safe
717
+		ParameterList thread_params;
718
+		#pragma omp critical
719
+		{
720
+			thread_params = ParameterList(params);
721
+		}
722
+
723
+		// update start jet
724
+		thread_params.sequence_start = thread_params.sequence_start + start_jet * Jets * steps * skip;
725
+
726
+		int start_format = (thread_params.file.find_last_of('/') + 1);
727
+		int end_format = thread_params.file.length() - start_format;
728
+
729
+		string sequence_path = thread_params.file.substr(0, start_format);
730
+
731
+		string format = "frame_%i.png";
732
+		string flow_format = thread_params.parameter<string>("flow_format", "frame_%i");
733
+
734
+		int len_format = flow_format.find_last_of('.');
735
+		flow_format = flow_format.substr(0, len_format);
736
+		format = thread_params.file.substr(start_format, end_format);
737
+		if (sequence_path[sequence_path.length() - 1] != '/')
738
+			sequence_path = sequence_path + "/";
739
+
740
+
741
+		char final_file[1024];
742
+		if (!sintel)
743
+			sprintf(final_file, (acc_folder.str() + "/" + flow_format + ".flo").c_str(), thread_params.sequence_start);
744
+		else
745
+			sprintf(final_file, (acc_folder.str() + "/s" + flow_format + ".flo").c_str(), thread_params.sequence_start, 0);
746
+
747
+		if (access(final_file, F_OK) != -1) {
748
+			cout << "Flow file " << final_file << " already exists!" << endl;
749
+			continue;
750
+		}
751
+		cout << "Flow file " << final_file << " does not exists!" << endl;
752
+
753
+		/*
754
+		 * ################### read in image sequence ###################
755
+		 */
756
+		vector<int> red_loc = thread_params.splitParameter<int>("raw_red_loc", "0,0");
757
+
758
+		// read in sequence, forward and backward flow, and segmentation
759
+		Mat *sequence = new Mat[Jets + 1];  						// high speed video sequence
760
+		Mat reference;
761
+		color_image_t **data = new color_image_t*[Jets + 1];  		// high speed video sequence
762
+		color_image_t **data_dx = new color_image_t*[Jets + 1];  	// derivative of high speed video sequence
763
+		color_image_t **data_dy = new color_image_t*[Jets + 1];  	// derivative of high speed video sequence
764
+		Mat* gt = NULL;
765
+		image_t*** gt_sfr = NULL;
766
+		Mat* gt_occlusions = NULL;
767
+		Mat* occlusions = new Mat[Jets];
768
+
769
+		Mat *forward_flow = new Mat[Jets];
770
+		Mat *backward_flow = new Mat[Jets];
771
+		/*
772
+		 * ################### read in image sequence ###################
773
+		 */
774
+		float norm = 1;
775
+		for (uint32_t f = 0; f < (Jets + 1); f++) {
776
+			char img_file[1024];
777
+			if (!sintel) {
778
+				sprintf(img_file, (sequence_path + format).c_str(), thread_params.sequence_start + f * steps * skip);
779
+			} else {
780
+				int sintel_frame = thread_params.sequence_start / 1000;
781
+				int hfr_frame = f * steps * skip + (thread_params.sequence_start % 1000);
782
+
783
+				while (hfr_frame < 0) {
784
+					sintel_frame--;
785
+					hfr_frame = 42 + hfr_frame;
786
+				}
787
+				while (hfr_frame > 41) {
788
+					sintel_frame++;
789
+					hfr_frame = hfr_frame - 42;
790
+				}
791
+
792
+				sprintf(img_file, (sequence_path + format).c_str(), sintel_frame, hfr_frame);
793
+			}
794
+			cout << "Reading " << img_file << "..." << endl;
795
+
796
+			sequence[f] = imread(string(img_file), CV_LOAD_IMAGE_UNCHANGED);         // load images
797
+
798
+			if (sequence[f].type() == 2 || sequence[f].type() == 18)
799
+				norm = 1.0f / 255;		// for 16 bit images
800
+
801
+			// convert to floating point
802
+			sequence[f].convertTo(sequence[f], CV_32FC(sequence[f].channels()));
803
+
804
+			/*
805
+			 * DEMOSAICING
806
+			 */
807
+			if (thread_params.exists("raw") && thread_params.parameter<bool>("raw")) {
808
+				Mat tmp = sequence[f].clone();
809
+				color_image_t* tmp_in = color_image_new(sequence[f].cols, sequence[f].rows);
810
+				color_image_t* tmp_out = color_image_new(sequence[f].cols, sequence[f].rows);
811
+
812
+				switch (thread_params.parameter<int>("raw_demosaicing", "0")) {
813
+					case 0: // use bilinear demosaicing
814
+						sequence[f] = Mat::zeros(tmp.rows, tmp.cols, CV_32FC3);
815
+						//						bayer2rgb(tmp, sequence[f], green_start, blue_start); 	// red green
816
+						bayer2rgbGR(tmp, sequence[f], red_loc[0], red_loc[1]); // red green
817
+						break;
818
+
819
+					case 1:	// use hamilton adams demosaicing
820
+						mat2colorImg<float>(sequence[f], tmp_in);
821
+
822
+						HADemosaicing(tmp_out->c1, tmp_in->c1, tmp_in->width, tmp_in->height, red_loc[0], red_loc[1]); // Hamilton-Adams implemented by Pascal Getreuer
823
+
824
+						sequence[f] = Mat::zeros(sequence[f].rows, sequence[f].cols, CV_32FC3);
825
+						colorImg2colorMat<Vec3f>(tmp_out, sequence[f]);
826
+						break;
827
+
828
+					case 2: // use opencv demosaicing
829
+						tmp.convertTo(tmp, CV_8UC1);
830
+						sequence[f] = Mat::zeros(tmp.rows, tmp.cols, CV_8UC3);
831
+
832
+						int code = CV_BayerBG2RGB;
833
+						if (red_loc[1] == 0) // y
834
+							if (red_loc[0] == 0) // x
835
+								code = CV_BayerBG2RGB;
836
+							else
837
+								code = CV_BayerGB2RGB;
838
+						else if (red_loc[0] == 0) // x
839
+							code = CV_BayerGR2RGB;
840
+						else
841
+							code = CV_BayerRG2RGB;
842
+
843
+						cv::cvtColor(tmp, sequence[f], code); // components from second row, second column !!!!!!!!!!!!!!!!!
844
+						sequence[f].convertTo(sequence[f], CV_32FC(sequence[f].channels()));
845
+						break;
846
+				}
847
+
848
+				color_image_delete(tmp_in);
849
+				color_image_delete(tmp_out);
850
+			} else {
851
+				// covert to RGB
852
+				cv::cvtColor(sequence[f], sequence[f], CV_BGR2RGB);
853
+			}
854
+
855
+			if (thread_params.parameter<bool>("grayscale", "0"))
856
+				cvtColor(sequence[f], sequence[f], CV_RGB2GRAY);
857
+
858
+			// use only a part of the images
859
+			if (thread_params.extent.x > 0 || thread_params.extent.y > 0) {
860
+				sequence[f] = sequence[f].rowRange(Range(thread_params.center.y - thread_params.extent.y / 2, thread_params.center.y + thread_params.extent.y / 2));
861
+				sequence[f] = sequence[f].colRange(Range(thread_params.center.x - thread_params.extent.x / 2, thread_params.center.x + thread_params.extent.x / 2));
862
+			}
863
+
864
+			// rescale image with gaussian blur to avoid anti-aliasing
865
+			double img_scale = thread_params.parameter<double>("scale", "1.0");
866
+			if (img_scale != 1) {
867
+				GaussianBlur(sequence[f], sequence[f], Size(0, 0), 1 / sqrt(2 * img_scale), 1 / sqrt(2 * img_scale), BORDER_REPLICATE);
868
+				resize(sequence[f], sequence[f], Size(0, 0), img_scale, img_scale, INTER_LINEAR);
869
+			}
870
+
871
+			// print to file
872
+			char file[1024];
873
+			sprintf(file, (acc_folder.str() + "sequence/frame_%i.png").c_str(), thread_params.sequence_start - steps * skip + f * skip);
874
+
875
+			Mat output_img;
876
+			if (thread_params.parameter<bool>("16bit", "0")) {
877
+				sequence[f].convertTo(output_img, CV_16UC(sequence[f].channels()));
878
+			} else {
879
+				sequence[f].convertTo(output_img, CV_8UC(sequence[f].channels()), norm);
880
+			}
881
+
882
+			if (thread_params.parameter<bool>("grayscale", "0"))
883
+				cv::cvtColor(output_img, output_img, CV_GRAY2BGR);	// OpenCV uses BGR
884
+			else
885
+				cv::cvtColor(output_img, output_img, CV_RGB2BGR);	// OpenCV uses BGR
886
+
887
+			if (thread_params.verbosity(WRITE_FILES)) {
888
+				imwrite(file, output_img, compression_params);
889
+			}
890
+
891
+			data[f] = color_image_new(sequence[f].cols, sequence[f].rows);
892
+			if (thread_params.parameter<bool>("grayscale", "0"))
893
+				mat2colorImg<float>(sequence[f], data[f]);
894
+			else
895
+				colorMat2colorImg<Vec3f>(sequence[f], data[f]);
896
+		}
897
+
898
+		// normalize data terms
899
+		normalize(data, Jets + 1, thread_params);
900
+
901
+		// compute derivatives
902
+		for (uint32_t f = 0; f < (uint32_t) (Jets + 1); f++) {
903
+			data_dx[f] = color_image_new(sequence[f].cols, sequence[f].rows);
904
+			data_dy[f] = color_image_new(sequence[f].cols, sequence[f].rows);
905
+
906
+			float deriv_filter[3] = { 0.0f, -8.0f / 12.0f, 1.0f / 12.0f };
907
+			convolution_t *deriv = convolution_new(2, deriv_filter, 0);
908
+			color_image_convolve_hv(data_dx[f], data[f], deriv, NULL);
909
+			color_image_convolve_hv(data_dy[f], data[f], NULL, deriv);
910
+			convolution_delete(deriv);
911
+		}
912
+
913
+		// store reference image
914
+		sequence[0].convertTo(reference, CV_8UC(sequence[0].channels()), norm);
915
+		double img_scale = 1.0f / (skip_pixel + 1);
916
+		if (img_scale != 1) {
917
+			GaussianBlur(reference, reference, Size(0, 0), 1 / sqrt(2 * img_scale), 1 / sqrt(2 * img_scale), BORDER_REPLICATE);
918
+			resize(reference, reference, Size(0, 0), img_scale, img_scale, INTER_LINEAR);
919
+		}
920
+
921
+		// reference image and edges for epic flow
922
+		float_image forward_edges;