4 // Code for computing optical flow between two images, and using it to interpolate
5 // in-between frames. The main user interface is the DISComputeFlow and Interpolate
6 // classes (also GrayscaleConversion can be useful).
17 // Predefined operating points from the paper.
18 struct OperatingPoint {
19 unsigned coarsest_level; // TODO: Adjust dynamically based on the resolution?
20 unsigned finest_level;
21 unsigned search_iterations; // Halved from the paper.
22 unsigned patch_size_pixels;
23 float patch_overlap_ratio;
24 bool variational_refinement;
26 // Not part of the original paper; used for interpolation.
27 // NOTE: Values much larger than 1.0 seems to trigger Haswell's “PMA stall”;
28 // the problem is not present on Broadwell and higher (there's a mitigation
29 // in the hardware, but Mesa doesn't enable it at the time of writing).
30 // Since we have hole filling, the holes from 1.0 are not critical,
31 // but larger values seem to do better than hole filling for large
32 // motion, blurs etc. since we have more candidates.
36 // Operating point 1 (600 Hz on CPU, excluding preprocessing).
37 static constexpr OperatingPoint operating_point1 = {
40 8, // Search iterations.
41 8, // Patch size (pixels).
42 0.30f, // Overlap ratio.
43 false, // Variational refinement.
44 1.0f // Splat size (pixels).
47 // Operating point 2 (300 Hz on CPU, excluding preprocessing).
48 static constexpr OperatingPoint operating_point2 = {
51 6, // Search iterations.
52 8, // Patch size (pixels).
53 0.40f, // Overlap ratio.
54 true, // Variational refinement.
55 1.0f // Splat size (pixels).
58 // Operating point 3 (10 Hz on CPU, excluding preprocessing).
59 // This is the only one that has been thorougly tested.
60 static constexpr OperatingPoint operating_point3 = {
63 8, // Search iterations.
64 12, // Patch size (pixels).
65 0.75f, // Overlap ratio.
66 true, // Variational refinement.
67 4.0f // Splat size (pixels).
70 // Operating point 4 (0.5 Hz on CPU, excluding preprocessing).
71 static constexpr OperatingPoint operating_point4 = {
74 128, // Search iterations.
75 12, // Patch size (pixels).
76 0.75f, // Overlap ratio.
77 true, // Variational refinement.
78 8.0f // Splat size (pixels).
81 int find_num_levels(int width, int height);
83 // A class that caches FBOs that render to a given set of textures.
84 // It never frees anything, so it is only suitable for rendering to
85 // the same (small) set of textures over and over again.
86 template<size_t num_elements>
87 class PersistentFBOSet {
89 void render_to(const std::array<GLuint, num_elements> &textures);
91 // Convenience wrappers.
92 void render_to(GLuint texture0) {
93 render_to({{texture0}});
96 void render_to(GLuint texture0, GLuint texture1) {
97 render_to({{texture0, texture1}});
100 void render_to(GLuint texture0, GLuint texture1, GLuint texture2) {
101 render_to({{texture0, texture1, texture2}});
104 void render_to(GLuint texture0, GLuint texture1, GLuint texture2, GLuint texture3) {
105 render_to({{texture0, texture1, texture2, texture3}});
109 // TODO: Delete these on destruction.
110 std::map<std::array<GLuint, num_elements>, GLuint> fbos;
114 // Same, but with a depth texture.
115 template<size_t num_elements>
116 class PersistentFBOSetWithDepth {
118 void render_to(GLuint depth_rb, const std::array<GLuint, num_elements> &textures);
120 // Convenience wrappers.
121 void render_to(GLuint depth_rb, GLuint texture0) {
122 render_to(depth_rb, {{texture0}});
125 void render_to(GLuint depth_rb, GLuint texture0, GLuint texture1) {
126 render_to(depth_rb, {{texture0, texture1}});
129 void render_to(GLuint depth_rb, GLuint texture0, GLuint texture1, GLuint texture2) {
130 render_to(depth_rb, {{texture0, texture1, texture2}});
133 void render_to(GLuint depth_rb, GLuint texture0, GLuint texture1, GLuint texture2, GLuint texture3) {
134 render_to(depth_rb, {{texture0, texture1, texture2, texture3}});
138 // TODO: Delete these on destruction.
139 std::map<std::pair<GLuint, std::array<GLuint, num_elements>>, GLuint> fbos;
142 // Convert RGB to grayscale, using Rec. 709 coefficients.
143 class GrayscaleConversion {
145 GrayscaleConversion();
146 void exec(GLint tex, GLint gray_tex, int width, int height, int num_layers);
149 PersistentFBOSet<1> fbos;
158 // Compute gradients in every point, used for the motion search.
159 // The DIS paper doesn't actually mention how these are computed,
160 // but seemingly, a 3x3 Sobel operator is used here (at least in
161 // later versions of the code), while a [1 -8 0 8 -1] kernel is
162 // used for all the derivatives in the variational refinement part
163 // (which borrows code from DeepFlow). This is inconsistent,
164 // but I guess we're better off with staying with the original
165 // decisions until we actually know having different ones would be better.
169 void exec(GLint tex_view, GLint grad_tex, int level_width, int level_height, int num_layers);
172 PersistentFBOSet<1> fbos;
175 GLuint sobel_program;
180 // Motion search to find the initial flow. See motion_search.frag for documentation.
183 MotionSearch(const OperatingPoint &op);
184 void exec(GLuint tex_view, GLuint grad_tex, GLuint flow_tex, GLuint flow_out_tex, int level_width, int level_height, int prev_level_width, int prev_level_height, int width_patches, int height_patches, int num_layers);
187 const OperatingPoint op;
188 PersistentFBOSet<1> fbos;
190 GLuint motion_vs_obj;
191 GLuint motion_fs_obj;
192 GLuint motion_search_program;
194 GLuint uniform_inv_image_size, uniform_inv_prev_level_size, uniform_out_flow_size;
195 GLuint uniform_image_tex, uniform_grad_tex, uniform_flow_tex;
196 GLuint uniform_patch_size, uniform_num_iterations;
199 // Do “densification”, ie., upsampling of the flow patches to the flow field
200 // (the same size as the image at this level). We draw one quad per patch
201 // over its entire covered area (using instancing in the vertex shader),
202 // and then weight the contributions in the pixel shader by post-warp difference.
203 // This is equation (3) in the paper.
205 // We accumulate the flow vectors in the R/G channels (for u/v) and the total
206 // weight in the B channel. Dividing R and G by B gives the normalized values.
209 Densify(const OperatingPoint &op);
210 void exec(GLuint tex_view, GLuint flow_tex, GLuint dense_flow_tex, int level_width, int level_height, int width_patches, int height_patches, int num_layers);
214 PersistentFBOSet<1> fbos;
216 GLuint densify_vs_obj;
217 GLuint densify_fs_obj;
218 GLuint densify_program;
220 GLuint uniform_patch_size;
221 GLuint uniform_image_tex, uniform_flow_tex;
224 // Warp I_1 to I_w, and then compute the mean (I) and difference (I_t) of
225 // I_0 and I_w. The prewarping is what enables us to solve the variational
226 // flow for du,dv instead of u,v.
228 // Also calculates the normalized flow, ie. divides by z (this is needed because
229 // Densify works by additive blending) and multiplies by the image size.
231 // See variational_refinement.txt for more information.
235 void exec(GLuint tex_view, GLuint flow_tex, GLuint normalized_flow_tex, GLuint I_tex, GLuint I_t_tex, int level_width, int level_height, int num_layers);
238 PersistentFBOSet<3> fbos;
240 GLuint prewarp_vs_obj;
241 GLuint prewarp_fs_obj;
242 GLuint prewarp_program;
244 GLuint uniform_image_tex, uniform_flow_tex;
247 // From I, calculate the partial derivatives I_x and I_y. We use a four-tap
248 // central difference filter, since apparently, that's tradition (I haven't
249 // measured quality versus a more normal 0.5 (I[x+1] - I[x-1]).)
250 // The coefficients come from
252 // https://en.wikipedia.org/wiki/Finite_difference_coefficient
254 // Also computes β_0, since it depends only on I_x and I_y.
258 void exec(GLuint input_tex, GLuint I_x_y_tex, GLuint beta_0_tex, int level_width, int level_height, int num_layers);
261 PersistentFBOSet<2> fbos;
263 GLuint derivatives_vs_obj;
264 GLuint derivatives_fs_obj;
265 GLuint derivatives_program;
270 // Calculate the diffusivity for each pixels, g(x,y). Smoothness (s) will
271 // be calculated in the shaders on-the-fly by sampling in-between two
272 // neighboring g(x,y) pixels, plus a border tweak to make sure we get
273 // zero smoothness at the border.
275 // See variational_refinement.txt for more information.
276 class ComputeDiffusivity {
278 ComputeDiffusivity();
279 void exec(GLuint flow_tex, GLuint diff_flow_tex, GLuint diffusivity_tex, int level_width, int level_height, bool zero_diff_flow, int num_layers);
282 PersistentFBOSet<1> fbos;
284 GLuint diffusivity_vs_obj;
285 GLuint diffusivity_fs_obj;
286 GLuint diffusivity_program;
288 GLuint uniform_flow_tex, uniform_diff_flow_tex;
289 GLuint uniform_alpha, uniform_zero_diff_flow;
292 // Set up the equations set (two equations in two unknowns, per pixel).
293 // We store five floats; the three non-redundant elements of the 2x2 matrix (A)
294 // as 32-bit floats, and the two elements on the right-hand side (b) as 16-bit
295 // floats. (Actually, we store the inverse of the diagonal elements, because
296 // we only ever need to divide by them.) This fits into four u32 values;
297 // R, G, B for the matrix (the last element is symmetric) and A for the two b values.
298 // All the values of the energy term (E_I, E_G, E_S), except the smoothness
299 // terms that depend on other pixels, are calculated in one pass.
301 // The equation set is split in two; one contains only the pixels needed for
302 // the red pass, and one only for the black pass (see sor.frag). This reduces
303 // the amount of data the SOR shader has to pull in, at the cost of some
304 // complexity when the equation texture ends up with half the size and we need
305 // to adjust texture coordinates. The contraction is done along the horizontal
306 // axis, so that on even rows (0, 2, 4, ...), the “red” texture will contain
307 // pixels 0, 2, 4, 6, etc., and on odd rows 1, 3, 5, etc..
309 // See variational_refinement.txt for more information about the actual
311 class SetupEquations {
314 void exec(GLuint I_x_y_tex, GLuint I_t_tex, GLuint diff_flow_tex, GLuint flow_tex, GLuint beta_0_tex, GLuint diffusivity_tex, GLuint equation_red_tex, GLuint equation_black_tex, int level_width, int level_height, bool zero_diff_flow, int num_layers);
317 PersistentFBOSet<2> fbos;
319 GLuint equations_vs_obj;
320 GLuint equations_fs_obj;
321 GLuint equations_program;
323 GLuint uniform_I_x_y_tex, uniform_I_t_tex;
324 GLuint uniform_diff_flow_tex, uniform_base_flow_tex;
325 GLuint uniform_beta_0_tex;
326 GLuint uniform_diffusivity_tex;
327 GLuint uniform_gamma, uniform_delta, uniform_zero_diff_flow;
330 // Actually solve the equation sets made by SetupEquations, by means of
331 // successive over-relaxation (SOR).
333 // See variational_refinement.txt for more information.
337 void exec(GLuint diff_flow_tex, GLuint equation_red_tex, GLuint equation_black_tex, GLuint diffusivity_tex, int level_width, int level_height, int num_iterations, bool zero_diff_flow, int num_layers, ScopedTimer *sor_timer);
340 PersistentFBOSet<1> fbos;
346 GLuint uniform_diff_flow_tex;
347 GLuint uniform_equation_red_tex, uniform_equation_black_tex;
348 GLuint uniform_diffusivity_tex;
349 GLuint uniform_phase, uniform_num_nonzero_phases;
352 // Simply add the differential flow found by the variational refinement to the base flow.
353 // The output is in base_flow_tex; we don't need to make a new texture.
357 void exec(GLuint base_flow_tex, GLuint diff_flow_tex, int level_width, int level_height, int num_layers);
360 PersistentFBOSet<1> fbos;
362 GLuint add_flow_vs_obj;
363 GLuint add_flow_fs_obj;
364 GLuint add_flow_program;
366 GLuint uniform_diff_flow_tex;
369 // Take a copy of the flow, bilinearly interpolated and scaled up.
373 void exec(GLuint in_tex, GLuint out_tex, int input_width, int input_height, int output_width, int output_height, int num_layers);
376 PersistentFBOSet<1> fbos;
378 GLuint resize_flow_vs_obj;
379 GLuint resize_flow_fs_obj;
380 GLuint resize_flow_program;
382 GLuint uniform_flow_tex;
383 GLuint uniform_scale_factor;
388 GLuint get_texture(GLenum format, GLuint width, GLuint height, GLuint num_layers = 0);
389 void release_texture(GLuint tex_num);
390 GLuint get_renderbuffer(GLenum format, GLuint width, GLuint height);
391 void release_renderbuffer(GLuint tex_num);
397 GLuint width, height, num_layers;
399 bool is_renderbuffer = false;
401 std::vector<Texture> textures;
404 class DISComputeFlow {
406 DISComputeFlow(int width, int height, const OperatingPoint &op);
412 enum ResizeStrategy {
414 RESIZE_FLOW_TO_FULL_SIZE
417 // The texture must have two layers (first and second frame).
418 // Returns a texture that must be released with release_texture()
420 GLuint exec(GLuint tex, FlowDirection flow_direction, ResizeStrategy resize_strategy);
422 void release_texture(GLuint tex) {
423 pool.release_texture(tex);
428 GLuint initial_flow_tex;
429 GLuint vertex_vbo, vao;
431 const OperatingPoint op;
433 // The various passes.
435 MotionSearch motion_search;
438 Derivatives derivatives;
439 ComputeDiffusivity compute_diffusivity;
440 SetupEquations setup_equations;
442 AddBaseFlow add_base_flow;
443 ResizeFlow resize_flow;
446 // Forward-warp the flow half-way (or rather, by alpha). A non-zero “splatting”
447 // radius fills most of the holes.
450 Splat(const OperatingPoint &op);
452 // alpha is the time of the interpolated frame (0..1).
453 void exec(GLuint image_tex, GLuint bidirectional_flow_tex, GLuint flow_tex, GLuint depth_rb, int width, int height, float alpha);
456 const OperatingPoint op;
457 PersistentFBOSetWithDepth<1> fbos;
461 GLuint splat_program;
463 GLuint uniform_splat_size, uniform_alpha;
464 GLuint uniform_image_tex, uniform_flow_tex;
465 GLuint uniform_inv_flow_size;
468 // Doing good and fast hole-filling on a GPU is nontrivial. We choose an option
469 // that's fairly simple (given that most holes are really small) and also hopefully
470 // cheap should the holes not be so small. Conceptually, we look for the first
471 // non-hole to the left of us (ie., shoot a ray until we hit something), then
472 // the first non-hole to the right of us, then up and down, and then average them
473 // all together. It's going to create “stars” if the holes are big, but OK, that's
476 // Our implementation here is efficient assuming that the hierarchical Z-buffer is
477 // on even for shaders that do discard (this typically kills early Z, but hopefully
478 // not hierarchical Z); we set up Z so that only holes are written to, which means
479 // that as soon as a hole is filled, the rasterizer should just skip it. Most of the
480 // fullscreen quads should just be discarded outright, really.
485 // Output will be in flow_tex, temp_tex[0, 1, 2], representing the filling
486 // from the down, left, right and up, respectively. Use HoleBlend to merge
488 void exec(GLuint flow_tex, GLuint depth_rb, GLuint temp_tex[3], int width, int height);
491 PersistentFBOSetWithDepth<1> fbos;
498 GLuint uniform_z, uniform_sample_offset;
501 // Blend the four directions from HoleFill into one pixel, so that single-pixel
502 // holes become the average of their four neighbors.
507 void exec(GLuint flow_tex, GLuint depth_rb, GLuint temp_tex[3], int width, int height);
510 PersistentFBOSetWithDepth<1> fbos;
514 GLuint blend_program;
516 GLuint uniform_left_tex, uniform_right_tex, uniform_up_tex, uniform_down_tex;
517 GLuint uniform_z, uniform_sample_offset;
523 void exec(GLuint image_tex, GLuint flow_tex, GLuint output_tex, int width, int height, float alpha);
526 PersistentFBOSet<1> fbos;
529 GLuint blend_program;
531 GLuint uniform_image_tex, uniform_flow_tex;
532 GLuint uniform_alpha, uniform_flow_consistency_tolerance;
537 Interpolate(int width, int height, const OperatingPoint &op);
539 // Returns a texture that must be released with release_texture()
540 // after use. image_tex must be a two-layer RGBA8 texture with mipmaps
541 // (unless flow_level == 0).
542 GLuint exec(GLuint image_tex, GLuint bidirectional_flow_tex, GLuint width, GLuint height, float alpha);
544 void release_texture(GLuint tex) {
545 pool.release_texture(tex);
549 int width, height, flow_level;
550 GLuint vertex_vbo, vao;
552 const OperatingPoint op;
556 HoleBlend hole_blend;
560 #endif // !defined(_FLOW_H)