git.sesse.net Git - nageru/blob - flow.cpp

   1 #define NO_SDL_GLEXT 1
   2
   3 #include <epoxy/gl.h>
   4
   5 #include <SDL2/SDL.h>
   6 #include <SDL2/SDL_error.h>
   7 #include <SDL2/SDL_events.h>
   8 #include <SDL2/SDL_image.h>
   9 #include <SDL2/SDL_keyboard.h>
  10 #include <SDL2/SDL_mouse.h>
  11 #include <SDL2/SDL_video.h>
  12
  13 #include <assert.h>
  14 #include <getopt.h>
  15 #include <stdio.h>
  16 #include <unistd.h>
  17
  18 #include "util.h"
  19
  20 #include <algorithm>
  21 #include <deque>
  22 #include <memory>
  23 #include <map>
  24 #include <stack>
  25 #include <vector>
  26
  27 #define BUFFER_OFFSET(i) ((char *)nullptr + (i))
  28
  29 using namespace std;
  30
  31 // Operating point 3 (10 Hz on CPU, excluding preprocessing).
  32 constexpr float patch_overlap_ratio = 0.75f;
  33 constexpr unsigned coarsest_level = 5;
  34 constexpr unsigned finest_level = 1;
  35 constexpr unsigned patch_size_pixels = 12;
  36
  37 // Weighting constants for the different parts of the variational refinement.
  38 // These don't correspond 1:1 to the values given in the DIS paper,
  39 // since we have different normalizations and ranges in some cases.
  40 // These are found through a simple grid search on some MPI-Sintel data,
  41 // although the error (EPE) seems to be fairly insensitive to the precise values.
  42 // Only the relative values matter, so we fix alpha (the smoothness constant)
  43 // at unity and tweak the others.
  44 float vr_alpha = 1.0f, vr_delta = 0.25f, vr_gamma = 0.25f;
  45
  46 bool enable_timing = true;
  47 bool enable_variational_refinement = true;  // Just for debugging.
  48
  49 // Some global OpenGL objects.
  50 // TODO: These should really be part of DISComputeFlow.
  51 GLuint nearest_sampler, linear_sampler, zero_border_sampler;
  52 GLuint vertex_vbo;
  53
  54 // Structures for asynchronous readback. We assume everything is the same size (and GL_RG16F).
  55 struct ReadInProgress {
  56         GLuint pbo;
  57         string filename0, filename1;
  58         string flow_filename, ppm_filename;  // Either may be empty for no write.
  59 };
  60 stack<GLuint> spare_pbos;
  61 deque<ReadInProgress> reads_in_progress;
  62
  63 int find_num_levels(int width, int height)
  64 {
  65         int levels = 1;
  66         for (int w = width, h = height; w > 1 || h > 1; ) {
  67                 w >>= 1;
  68                 h >>= 1;
  69                 ++levels;
  70         }
  71         return levels;
  72 }
  73
  74 string read_file(const string &filename)
  75 {
  76         FILE *fp = fopen(filename.c_str(), "r");
  77         if (fp == nullptr) {
  78                 perror(filename.c_str());
  79                 exit(1);
  80         }
  81
  82         int ret = fseek(fp, 0, SEEK_END);
  83         if (ret == -1) {
  84                 perror("fseek(SEEK_END)");
  85                 exit(1);
  86         }
  87
  88         int size = ftell(fp);
  89
  90         ret = fseek(fp, 0, SEEK_SET);
  91         if (ret == -1) {
  92                 perror("fseek(SEEK_SET)");
  93                 exit(1);
  94         }
  95
  96         string str;
  97         str.resize(size);
  98         ret = fread(&str[0], size, 1, fp);
  99         if (ret == -1) {
 100                 perror("fread");
 101                 exit(1);
 102         }
 103         if (ret == 0) {
 104                 fprintf(stderr, "Short read when trying to read %d bytes from %s\n",
 105                                 size, filename.c_str());
 106                 exit(1);
 107         }
 108         fclose(fp);
 109
 110         return str;
 111 }
 112
 113
 114 GLuint compile_shader(const string &shader_src, GLenum type)
 115 {
 116         GLuint obj = glCreateShader(type);
 117         const GLchar* source[] = { shader_src.data() };
 118         const GLint length[] = { (GLint)shader_src.size() };
 119         glShaderSource(obj, 1, source, length);
 120         glCompileShader(obj);
 121
 122         GLchar info_log[4096];
 123         GLsizei log_length = sizeof(info_log) - 1;
 124         glGetShaderInfoLog(obj, log_length, &log_length, info_log);
 125         info_log[log_length] = 0;
 126         if (strlen(info_log) > 0) {
 127                 fprintf(stderr, "Shader compile log: %s\n", info_log);
 128         }
 129
 130         GLint status;
 131         glGetShaderiv(obj, GL_COMPILE_STATUS, &status);
 132         if (status == GL_FALSE) {
 133                 // Add some line numbers to easier identify compile errors.
 134                 string src_with_lines = "/*   1 */ ";
 135                 size_t lineno = 1;
 136                 for (char ch : shader_src) {
 137                         src_with_lines.push_back(ch);
 138                         if (ch == '\n') {
 139                                 char buf[32];
 140                                 snprintf(buf, sizeof(buf), "/* %3zu */ ", ++lineno);
 141                                 src_with_lines += buf;
 142                         }
 143                 }
 144
 145                 fprintf(stderr, "Failed to compile shader:\n%s\n", src_with_lines.c_str());
 146                 exit(1);
 147         }
 148
 149         return obj;
 150 }
 151
 152 GLuint load_texture(const char *filename, unsigned *width_ret, unsigned *height_ret)
 153 {
 154         SDL_Surface *surf = IMG_Load(filename);
 155         if (surf == nullptr) {
 156                 fprintf(stderr, "IMG_Load(%s): %s\n", filename, IMG_GetError());
 157                 exit(1);
 158         }
 159
 160         // For whatever reason, SDL doesn't support converting to YUV surfaces
 161         // nor grayscale, so we'll do it ourselves.
 162         SDL_Surface *rgb_surf = SDL_ConvertSurfaceFormat(surf, SDL_PIXELFORMAT_RGBA32, /*flags=*/0);
 163         if (rgb_surf == nullptr) {
 164                 fprintf(stderr, "SDL_ConvertSurfaceFormat(%s): %s\n", filename, SDL_GetError());
 165                 exit(1);
 166         }
 167
 168         SDL_FreeSurface(surf);
 169
 170         unsigned width = rgb_surf->w, height = rgb_surf->h;
 171         const uint8_t *sptr = (uint8_t *)rgb_surf->pixels;
 172         unique_ptr<uint8_t[]> pix(new uint8_t[width * height * 4]);
 173
 174         // Extract the Y component, and convert to bottom-left origin.
 175         for (unsigned y = 0; y < height; ++y) {
 176                 unsigned y2 = height - 1 - y;
 177                 memcpy(pix.get() + y * width * 4, sptr + y2 * rgb_surf->pitch, width * 4);
 178         }
 179         SDL_FreeSurface(rgb_surf);
 180
 181         GLuint tex;
 182         glCreateTextures(GL_TEXTURE_2D, 1, &tex);
 183         glTextureStorage2D(tex, 1, GL_RGBA8, width, height);
 184         glTextureSubImage2D(tex, 0, 0, 0, width, height, GL_RGBA, GL_UNSIGNED_BYTE, pix.get());
 185
 186         *width_ret = width;
 187         *height_ret = height;
 188
 189         return tex;
 190 }
 191
 192 GLuint link_program(GLuint vs_obj, GLuint fs_obj)
 193 {
 194         GLuint program = glCreateProgram();
 195         glAttachShader(program, vs_obj);
 196         glAttachShader(program, fs_obj);
 197         glLinkProgram(program);
 198         GLint success;
 199         glGetProgramiv(program, GL_LINK_STATUS, &success);
 200         if (success == GL_FALSE) {
 201                 GLchar error_log[1024] = {0};
 202                 glGetProgramInfoLog(program, 1024, nullptr, error_log);
 203                 fprintf(stderr, "Error linking program: %s\n", error_log);
 204                 exit(1);
 205         }
 206         return program;
 207 }
 208
 209 GLuint generate_vbo(GLint size, GLsizeiptr data_size, const GLvoid *data)
 210 {
 211         GLuint vbo;
 212         glCreateBuffers(1, &vbo);
 213         glBufferData(GL_ARRAY_BUFFER, data_size, data, GL_STATIC_DRAW);
 214         glNamedBufferData(vbo, data_size, data, GL_STATIC_DRAW);
 215         return vbo;
 216 }
 217
 218 GLuint fill_vertex_attribute(GLuint vao, GLuint glsl_program_num, const string &attribute_name, GLint size, GLenum type, GLsizeiptr data_size, const GLvoid *data)
 219 {
 220         int attrib = glGetAttribLocation(glsl_program_num, attribute_name.c_str());
 221         if (attrib == -1) {
 222                 return -1;
 223         }
 224
 225         GLuint vbo = generate_vbo(size, data_size, data);
 226
 227         glBindBuffer(GL_ARRAY_BUFFER, vbo);
 228         glEnableVertexArrayAttrib(vao, attrib);
 229         glVertexAttribPointer(attrib, size, type, GL_FALSE, 0, BUFFER_OFFSET(0));
 230         glBindBuffer(GL_ARRAY_BUFFER, 0);
 231
 232         return vbo;
 233 }
 234
 235 void bind_sampler(GLuint program, GLint location, GLuint texture_unit, GLuint tex, GLuint sampler)
 236 {
 237         if (location == -1) {
 238                 return;
 239         }
 240
 241         glBindTextureUnit(texture_unit, tex);
 242         glBindSampler(texture_unit, sampler);
 243         glProgramUniform1i(program, location, texture_unit);
 244 }
 245
 246 // A class that caches FBOs that render to a given set of textures.
 247 // It never frees anything, so it is only suitable for rendering to
 248 // the same (small) set of textures over and over again.
 249 template<size_t num_elements>
 250 class PersistentFBOSet {
 251 public:
 252         void render_to(const array<GLuint, num_elements> &textures);
 253
 254         // Convenience wrappers.
 255         void render_to(GLuint texture0, enable_if<num_elements == 1> * = nullptr) {
 256                 render_to({{texture0}});
 257         }
 258
 259         void render_to(GLuint texture0, GLuint texture1, enable_if<num_elements == 2> * = nullptr) {
 260                 render_to({{texture0, texture1}});
 261         }
 262
 263         void render_to(GLuint texture0, GLuint texture1, GLuint texture2, enable_if<num_elements == 3> * = nullptr) {
 264                 render_to({{texture0, texture1, texture2}});
 265         }
 266
 267         void render_to(GLuint texture0, GLuint texture1, GLuint texture2, GLuint texture3, enable_if<num_elements == 4> * = nullptr) {
 268                 render_to({{texture0, texture1, texture2, texture3}});
 269         }
 270
 271 private:
 272         // TODO: Delete these on destruction.
 273         map<array<GLuint, num_elements>, GLuint> fbos;
 274 };
 275
 276 template<size_t num_elements>
 277 void PersistentFBOSet<num_elements>::render_to(const array<GLuint, num_elements> &textures)
 278 {
 279         auto it = fbos.find(textures);
 280         if (it != fbos.end()) {
 281                 glBindFramebuffer(GL_FRAMEBUFFER, it->second);
 282                 return;
 283         }
 284
 285         GLuint fbo;
 286         glCreateFramebuffers(1, &fbo);
 287         GLenum bufs[num_elements];
 288         for (size_t i = 0; i < num_elements; ++i) {
 289                 glNamedFramebufferTexture(fbo, GL_COLOR_ATTACHMENT0 + i, textures[i], 0);
 290                 bufs[i] = GL_COLOR_ATTACHMENT0 + i;
 291         }
 292         glNamedFramebufferDrawBuffers(fbo, num_elements, bufs);
 293
 294         fbos[textures] = fbo;
 295         glBindFramebuffer(GL_FRAMEBUFFER, fbo);
 296 }
 297
 298 // Convert RGB to grayscale, using Rec. 709 coefficients.
 299 class GrayscaleConversion {
 300 public:
 301         GrayscaleConversion();
 302         void exec(GLint tex, GLint gray_tex, int width, int height);
 303
 304 private:
 305         PersistentFBOSet<1> fbos;
 306         GLuint gray_vs_obj;
 307         GLuint gray_fs_obj;
 308         GLuint gray_program;
 309         GLuint gray_vao;
 310
 311         GLuint uniform_tex;
 312 };
 313
 314 GrayscaleConversion::GrayscaleConversion()
 315 {
 316         gray_vs_obj = compile_shader(read_file("vs.vert"), GL_VERTEX_SHADER);
 317         gray_fs_obj = compile_shader(read_file("gray.frag"), GL_FRAGMENT_SHADER);
 318         gray_program = link_program(gray_vs_obj, gray_fs_obj);
 319
 320         // Set up the VAO containing all the required position/texcoord data.
 321         glCreateVertexArrays(1, &gray_vao);
 322         glBindVertexArray(gray_vao);
 323
 324         GLint position_attrib = glGetAttribLocation(gray_program, "position");
 325         glEnableVertexArrayAttrib(gray_vao, position_attrib);
 326         glVertexAttribPointer(position_attrib, 2, GL_FLOAT, GL_FALSE, 0, BUFFER_OFFSET(0));
 327
 328         uniform_tex = glGetUniformLocation(gray_program, "tex");
 329 }
 330
 331 void GrayscaleConversion::exec(GLint tex, GLint gray_tex, int width, int height)
 332 {
 333         glUseProgram(gray_program);
 334         bind_sampler(gray_program, uniform_tex, 0, tex, nearest_sampler);
 335
 336         glViewport(0, 0, width, height);
 337         fbos.render_to(gray_tex);
 338         glBindVertexArray(gray_vao);
 339         glUseProgram(gray_program);
 340         glDisable(GL_BLEND);
 341         glDrawArrays(GL_TRIANGLE_STRIP, 0, 4);
 342 }
 343
 344 // Compute gradients in every point, used for the motion search.
 345 // The DIS paper doesn't actually mention how these are computed,
 346 // but seemingly, a 3x3 Sobel operator is used here (at least in
 347 // later versions of the code), while a [1 -8 0 8 -1] kernel is
 348 // used for all the derivatives in the variational refinement part
 349 // (which borrows code from DeepFlow). This is inconsistent,
 350 // but I guess we're better off with staying with the original
 351 // decisions until we actually know having different ones would be better.
 352 class Sobel {
 353 public:
 354         Sobel();
 355         void exec(GLint tex0_view, GLint grad0_tex, int level_width, int level_height);
 356
 357 private:
 358         PersistentFBOSet<1> fbos;
 359         GLuint sobel_vs_obj;
 360         GLuint sobel_fs_obj;
 361         GLuint sobel_program;
 362         GLuint sobel_vao;
 363
 364         GLuint uniform_tex;
 365 };
 366
 367 Sobel::Sobel()
 368 {
 369         sobel_vs_obj = compile_shader(read_file("vs.vert"), GL_VERTEX_SHADER);
 370         sobel_fs_obj = compile_shader(read_file("sobel.frag"), GL_FRAGMENT_SHADER);
 371         sobel_program = link_program(sobel_vs_obj, sobel_fs_obj);
 372
 373         // Set up the VAO containing all the required position/texcoord data.
 374         glCreateVertexArrays(1, &sobel_vao);
 375         glBindVertexArray(sobel_vao);
 376
 377         GLint position_attrib = glGetAttribLocation(sobel_program, "position");
 378         glEnableVertexArrayAttrib(sobel_vao, position_attrib);
 379         glVertexAttribPointer(position_attrib, 2, GL_FLOAT, GL_FALSE, 0, BUFFER_OFFSET(0));
 380
 381         uniform_tex = glGetUniformLocation(sobel_program, "tex");
 382 }
 383
 384 void Sobel::exec(GLint tex0_view, GLint grad0_tex, int level_width, int level_height)
 385 {
 386         glUseProgram(sobel_program);
 387         bind_sampler(sobel_program, uniform_tex, 0, tex0_view, nearest_sampler);
 388
 389         glViewport(0, 0, level_width, level_height);
 390         fbos.render_to(grad0_tex);
 391         glBindVertexArray(sobel_vao);
 392         glUseProgram(sobel_program);
 393         glDisable(GL_BLEND);
 394         glDrawArrays(GL_TRIANGLE_STRIP, 0, 4);
 395 }
 396
 397 // Motion search to find the initial flow. See motion_search.frag for documentation.
 398 class MotionSearch {
 399 public:
 400         MotionSearch();
 401         void exec(GLuint tex0_view, GLuint tex1_view, GLuint grad0_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);
 402
 403 private:
 404         PersistentFBOSet<1> fbos;
 405
 406         GLuint motion_vs_obj;
 407         GLuint motion_fs_obj;
 408         GLuint motion_search_program;
 409         GLuint motion_search_vao;
 410
 411         GLuint uniform_inv_image_size, uniform_inv_prev_level_size;
 412         GLuint uniform_image0_tex, uniform_image1_tex, uniform_grad0_tex, uniform_flow_tex;
 413 };
 414
 415 MotionSearch::MotionSearch()
 416 {
 417         motion_vs_obj = compile_shader(read_file("motion_search.vert"), GL_VERTEX_SHADER);
 418         motion_fs_obj = compile_shader(read_file("motion_search.frag"), GL_FRAGMENT_SHADER);
 419         motion_search_program = link_program(motion_vs_obj, motion_fs_obj);
 420
 421         // Set up the VAO containing all the required position/texcoord data.
 422         glCreateVertexArrays(1, &motion_search_vao);
 423         glBindVertexArray(motion_search_vao);
 424         glBindBuffer(GL_ARRAY_BUFFER, vertex_vbo);
 425
 426         GLint position_attrib = glGetAttribLocation(motion_search_program, "position");
 427         glEnableVertexArrayAttrib(motion_search_vao, position_attrib);
 428         glVertexAttribPointer(position_attrib, 2, GL_FLOAT, GL_FALSE, 0, BUFFER_OFFSET(0));
 429
 430         uniform_inv_image_size = glGetUniformLocation(motion_search_program, "inv_image_size");
 431         uniform_inv_prev_level_size = glGetUniformLocation(motion_search_program, "inv_prev_level_size");
 432         uniform_image0_tex = glGetUniformLocation(motion_search_program, "image0_tex");
 433         uniform_image1_tex = glGetUniformLocation(motion_search_program, "image1_tex");
 434         uniform_grad0_tex = glGetUniformLocation(motion_search_program, "grad0_tex");
 435         uniform_flow_tex = glGetUniformLocation(motion_search_program, "flow_tex");
 436 }
 437
 438 void MotionSearch::exec(GLuint tex0_view, GLuint tex1_view, GLuint grad0_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)
 439 {
 440         glUseProgram(motion_search_program);
 441
 442         bind_sampler(motion_search_program, uniform_image0_tex, 0, tex0_view, nearest_sampler);
 443         bind_sampler(motion_search_program, uniform_image1_tex, 1, tex1_view, linear_sampler);
 444         bind_sampler(motion_search_program, uniform_grad0_tex, 2, grad0_tex, zero_border_sampler);
 445         bind_sampler(motion_search_program, uniform_flow_tex, 3, flow_tex, linear_sampler);
 446
 447         glProgramUniform2f(motion_search_program, uniform_inv_image_size, 1.0f / level_width, 1.0f / level_height);
 448         glProgramUniform2f(motion_search_program, uniform_inv_prev_level_size, 1.0f / prev_level_width, 1.0f / prev_level_height);
 449
 450         glViewport(0, 0, width_patches, height_patches);
 451         fbos.render_to(flow_out_tex);
 452         glBindVertexArray(motion_search_vao);
 453         glUseProgram(motion_search_program);
 454         glDrawArrays(GL_TRIANGLE_STRIP, 0, 4);
 455 }
 456
 457 // Do “densification”, ie., upsampling of the flow patches to the flow field
 458 // (the same size as the image at this level). We draw one quad per patch
 459 // over its entire covered area (using instancing in the vertex shader),
 460 // and then weight the contributions in the pixel shader by post-warp difference.
 461 // This is equation (3) in the paper.
 462 //
 463 // We accumulate the flow vectors in the R/G channels (for u/v) and the total
 464 // weight in the B channel. Dividing R and G by B gives the normalized values.
 465 class Densify {
 466 public:
 467         Densify();
 468         void exec(GLuint tex0_view, GLuint tex1_view, GLuint flow_tex, GLuint dense_flow_tex, int level_width, int level_height, int width_patches, int height_patches);
 469
 470 private:
 471         PersistentFBOSet<1> fbos;
 472
 473         GLuint densify_vs_obj;
 474         GLuint densify_fs_obj;
 475         GLuint densify_program;
 476         GLuint densify_vao;
 477
 478         GLuint uniform_patch_size, uniform_patch_spacing;
 479         GLuint uniform_image0_tex, uniform_image1_tex, uniform_flow_tex;
 480 };
 481
 482 Densify::Densify()
 483 {
 484         densify_vs_obj = compile_shader(read_file("densify.vert"), GL_VERTEX_SHADER);
 485         densify_fs_obj = compile_shader(read_file("densify.frag"), GL_FRAGMENT_SHADER);
 486         densify_program = link_program(densify_vs_obj, densify_fs_obj);
 487
 488         // Set up the VAO containing all the required position/texcoord data.
 489         glCreateVertexArrays(1, &densify_vao);
 490         glBindVertexArray(densify_vao);
 491         glBindBuffer(GL_ARRAY_BUFFER, vertex_vbo);
 492
 493         GLint position_attrib = glGetAttribLocation(densify_program, "position");
 494         glEnableVertexArrayAttrib(densify_vao, position_attrib);
 495         glVertexAttribPointer(position_attrib, 2, GL_FLOAT, GL_FALSE, 0, BUFFER_OFFSET(0));
 496
 497         uniform_patch_size = glGetUniformLocation(densify_program, "patch_size");
 498         uniform_patch_spacing = glGetUniformLocation(densify_program, "patch_spacing");
 499         uniform_image0_tex = glGetUniformLocation(densify_program, "image0_tex");
 500         uniform_image1_tex = glGetUniformLocation(densify_program, "image1_tex");
 501         uniform_flow_tex = glGetUniformLocation(densify_program, "flow_tex");
 502 }
 503
 504 void Densify::exec(GLuint tex0_view, GLuint tex1_view, GLuint flow_tex, GLuint dense_flow_tex, int level_width, int level_height, int width_patches, int height_patches)
 505 {
 506         glUseProgram(densify_program);
 507
 508         bind_sampler(densify_program, uniform_image0_tex, 0, tex0_view, nearest_sampler);
 509         bind_sampler(densify_program, uniform_image1_tex, 1, tex1_view, linear_sampler);
 510         bind_sampler(densify_program, uniform_flow_tex, 2, flow_tex, nearest_sampler);
 511
 512         glProgramUniform2f(densify_program, uniform_patch_size,
 513                 float(patch_size_pixels) / level_width,
 514                 float(patch_size_pixels) / level_height);
 515
 516         float patch_spacing_x = float(level_width - patch_size_pixels) / (width_patches - 1);
 517         float patch_spacing_y = float(level_height - patch_size_pixels) / (height_patches - 1);
 518         if (width_patches == 1) patch_spacing_x = 0.0f;  // Avoid infinities.
 519         if (height_patches == 1) patch_spacing_y = 0.0f;
 520         glProgramUniform2f(densify_program, uniform_patch_spacing,
 521                 patch_spacing_x / level_width,
 522                 patch_spacing_y / level_height);
 523
 524         glViewport(0, 0, level_width, level_height);
 525         glEnable(GL_BLEND);
 526         glBlendFunc(GL_ONE, GL_ONE);
 527         glBindVertexArray(densify_vao);
 528         fbos.render_to(dense_flow_tex);
 529         glDrawArraysInstanced(GL_TRIANGLE_STRIP, 0, 4, width_patches * height_patches);
 530 }
 531
 532 // Warp I_1 to I_w, and then compute the mean (I) and difference (I_t) of
 533 // I_0 and I_w. The prewarping is what enables us to solve the variational
 534 // flow for du,dv instead of u,v.
 535 //
 536 // Also calculates the normalized flow, ie. divides by z (this is needed because
 537 // Densify works by additive blending) and multiplies by the image size.
 538 //
 539 // See variational_refinement.txt for more information.
 540 class Prewarp {
 541 public:
 542         Prewarp();
 543         void exec(GLuint tex0_view, GLuint tex1_view, GLuint flow_tex, GLuint normalized_flow_tex, GLuint I_tex, GLuint I_t_tex, int level_width, int level_height);
 544
 545 private:
 546         PersistentFBOSet<3> fbos;
 547
 548         GLuint prewarp_vs_obj;
 549         GLuint prewarp_fs_obj;
 550         GLuint prewarp_program;
 551         GLuint prewarp_vao;
 552
 553         GLuint uniform_image0_tex, uniform_image1_tex, uniform_flow_tex;
 554 };
 555
 556 Prewarp::Prewarp()
 557 {
 558         prewarp_vs_obj = compile_shader(read_file("vs.vert"), GL_VERTEX_SHADER);
 559         prewarp_fs_obj = compile_shader(read_file("prewarp.frag"), GL_FRAGMENT_SHADER);
 560         prewarp_program = link_program(prewarp_vs_obj, prewarp_fs_obj);
 561
 562         // Set up the VAO containing all the required position/texcoord data.
 563         glCreateVertexArrays(1, &prewarp_vao);
 564         glBindVertexArray(prewarp_vao);
 565         glBindBuffer(GL_ARRAY_BUFFER, vertex_vbo);
 566
 567         GLint position_attrib = glGetAttribLocation(prewarp_program, "position");
 568         glEnableVertexArrayAttrib(prewarp_vao, position_attrib);
 569         glVertexAttribPointer(position_attrib, 2, GL_FLOAT, GL_FALSE, 0, BUFFER_OFFSET(0));
 570
 571         uniform_image0_tex = glGetUniformLocation(prewarp_program, "image0_tex");
 572         uniform_image1_tex = glGetUniformLocation(prewarp_program, "image1_tex");
 573         uniform_flow_tex = glGetUniformLocation(prewarp_program, "flow_tex");
 574 }
 575
 576 void Prewarp::exec(GLuint tex0_view, GLuint tex1_view, GLuint flow_tex, GLuint I_tex, GLuint I_t_tex, GLuint normalized_flow_tex, int level_width, int level_height)
 577 {
 578         glUseProgram(prewarp_program);
 579
 580         bind_sampler(prewarp_program, uniform_image0_tex, 0, tex0_view, nearest_sampler);
 581         bind_sampler(prewarp_program, uniform_image1_tex, 1, tex1_view, linear_sampler);
 582         bind_sampler(prewarp_program, uniform_flow_tex, 2, flow_tex, nearest_sampler);
 583
 584         glViewport(0, 0, level_width, level_height);
 585         glDisable(GL_BLEND);
 586         glBindVertexArray(prewarp_vao);
 587         fbos.render_to(I_tex, I_t_tex, normalized_flow_tex);
 588         glDrawArrays(GL_TRIANGLE_STRIP, 0, 4);
 589 }
 590
 591 // From I, calculate the partial derivatives I_x and I_y. We use a four-tap
 592 // central difference filter, since apparently, that's tradition (I haven't
 593 // measured quality versus a more normal 0.5 (I[x+1] - I[x-1]).)
 594 // The coefficients come from
 595 //
 596 //   https://en.wikipedia.org/wiki/Finite_difference_coefficient
 597 //
 598 // Also computes β_0, since it depends only on I_x and I_y.
 599 class Derivatives {
 600 public:
 601         Derivatives();
 602         void exec(GLuint input_tex, GLuint I_x_y_tex, GLuint beta_0_tex, int level_width, int level_height);
 603
 604 private:
 605         PersistentFBOSet<2> fbos;
 606
 607         GLuint derivatives_vs_obj;
 608         GLuint derivatives_fs_obj;
 609         GLuint derivatives_program;
 610         GLuint derivatives_vao;
 611
 612         GLuint uniform_tex;
 613 };
 614
 615 Derivatives::Derivatives()
 616 {
 617         derivatives_vs_obj = compile_shader(read_file("vs.vert"), GL_VERTEX_SHADER);
 618         derivatives_fs_obj = compile_shader(read_file("derivatives.frag"), GL_FRAGMENT_SHADER);
 619         derivatives_program = link_program(derivatives_vs_obj, derivatives_fs_obj);
 620
 621         // Set up the VAO containing all the required position/texcoord data.
 622         glCreateVertexArrays(1, &derivatives_vao);
 623         glBindVertexArray(derivatives_vao);
 624         glBindBuffer(GL_ARRAY_BUFFER, vertex_vbo);
 625
 626         GLint position_attrib = glGetAttribLocation(derivatives_program, "position");
 627         glEnableVertexArrayAttrib(derivatives_vao, position_attrib);
 628         glVertexAttribPointer(position_attrib, 2, GL_FLOAT, GL_FALSE, 0, BUFFER_OFFSET(0));
 629
 630         uniform_tex = glGetUniformLocation(derivatives_program, "tex");
 631 }
 632
 633 void Derivatives::exec(GLuint input_tex, GLuint I_x_y_tex, GLuint beta_0_tex, int level_width, int level_height)
 634 {
 635         glUseProgram(derivatives_program);
 636
 637         bind_sampler(derivatives_program, uniform_tex, 0, input_tex, nearest_sampler);
 638
 639         glViewport(0, 0, level_width, level_height);
 640         glDisable(GL_BLEND);
 641         glBindVertexArray(derivatives_vao);
 642         fbos.render_to(I_x_y_tex, beta_0_tex);
 643         glDrawArrays(GL_TRIANGLE_STRIP, 0, 4);
 644 }
 645
 646 // Calculate the smoothness constraints between neighboring pixels;
 647 // s_x(x,y) stores smoothness between pixel (x,y) and (x+1,y),
 648 // and s_y(x,y) stores between (x,y) and (x,y+1). We'll sample with
 649 // border color (0,0) later, so that there's zero diffusion out of
 650 // the border.
 651 //
 652 // See variational_refinement.txt for more information.
 653 class ComputeSmoothness {
 654 public:
 655         ComputeSmoothness();
 656         void exec(GLuint flow_tex, GLuint diff_flow_tex, GLuint smoothness_x_tex, GLuint smoothness_y_tex, int level_width, int level_height);
 657
 658 private:
 659         PersistentFBOSet<2> fbos;
 660
 661         GLuint smoothness_vs_obj;
 662         GLuint smoothness_fs_obj;
 663         GLuint smoothness_program;
 664         GLuint smoothness_vao;
 665
 666         GLuint uniform_flow_tex, uniform_diff_flow_tex;
 667         GLuint uniform_alpha;
 668 };
 669
 670 ComputeSmoothness::ComputeSmoothness()
 671 {
 672         smoothness_vs_obj = compile_shader(read_file("vs.vert"), GL_VERTEX_SHADER);
 673         smoothness_fs_obj = compile_shader(read_file("smoothness.frag"), GL_FRAGMENT_SHADER);
 674         smoothness_program = link_program(smoothness_vs_obj, smoothness_fs_obj);
 675
 676         // Set up the VAO containing all the required position/texcoord data.
 677         glCreateVertexArrays(1, &smoothness_vao);
 678         glBindVertexArray(smoothness_vao);
 679         glBindBuffer(GL_ARRAY_BUFFER, vertex_vbo);
 680
 681         GLint position_attrib = glGetAttribLocation(smoothness_program, "position");
 682         glEnableVertexArrayAttrib(smoothness_vao, position_attrib);
 683         glVertexAttribPointer(position_attrib, 2, GL_FLOAT, GL_FALSE, 0, BUFFER_OFFSET(0));
 684
 685         uniform_flow_tex = glGetUniformLocation(smoothness_program, "flow_tex");
 686         uniform_diff_flow_tex = glGetUniformLocation(smoothness_program, "diff_flow_tex");
 687         uniform_alpha = glGetUniformLocation(smoothness_program, "alpha");
 688 }
 689
 690 void ComputeSmoothness::exec(GLuint flow_tex, GLuint diff_flow_tex, GLuint smoothness_x_tex, GLuint smoothness_y_tex, int level_width, int level_height)
 691 {
 692         glUseProgram(smoothness_program);
 693
 694         bind_sampler(smoothness_program, uniform_flow_tex, 0, flow_tex, nearest_sampler);
 695         bind_sampler(smoothness_program, uniform_diff_flow_tex, 1, diff_flow_tex, nearest_sampler);
 696         glProgramUniform1f(smoothness_program, uniform_alpha, vr_alpha);
 697
 698         glViewport(0, 0, level_width, level_height);
 699
 700         glDisable(GL_BLEND);
 701         glBindVertexArray(smoothness_vao);
 702         fbos.render_to(smoothness_x_tex, smoothness_y_tex);
 703         glDrawArrays(GL_TRIANGLE_STRIP, 0, 4);
 704
 705         // Make sure the smoothness on the right and upper borders is zero.
 706         // We could have done this by making (W-1)xH and Wx(H-1) textures instead
 707         // (we're sampling smoothness with all-zero border color), but we'd
 708         // have to adjust the sampling coordinates, which is annoying.
 709         glClearTexSubImage(smoothness_x_tex, 0,  level_width - 1, 0, 0,   1, level_height, 1,  GL_RED, GL_FLOAT, nullptr);
 710         glClearTexSubImage(smoothness_y_tex, 0,  0, level_height - 1, 0,  level_width, 1, 1,   GL_RED, GL_FLOAT, nullptr);
 711 }
 712
 713 // Set up the equations set (two equations in two unknowns, per pixel).
 714 // We store five floats; the three non-redundant elements of the 2x2 matrix (A)
 715 // as 32-bit floats, and the two elements on the right-hand side (b) as 16-bit
 716 // floats. (Actually, we store the inverse of the diagonal elements, because
 717 // we only ever need to divide by them.) This fits into four u32 values;
 718 // R, G, B for the matrix (the last element is symmetric) and A for the two b values.
 719 // All the values of the energy term (E_I, E_G, E_S), except the smoothness
 720 // terms that depend on other pixels, are calculated in one pass.
 721 //
 722 // See variational_refinement.txt for more information.
 723 class SetupEquations {
 724 public:
 725         SetupEquations();
 726         void exec(GLuint I_x_y_tex, GLuint I_t_tex, GLuint diff_flow_tex, GLuint flow_tex, GLuint beta_0_tex, GLuint smoothness_x_tex, GLuint smoothness_y_tex, GLuint equation_tex, int level_width, int level_height);
 727
 728 private:
 729         PersistentFBOSet<1> fbos;
 730
 731         GLuint equations_vs_obj;
 732         GLuint equations_fs_obj;
 733         GLuint equations_program;
 734         GLuint equations_vao;
 735
 736         GLuint uniform_I_x_y_tex, uniform_I_t_tex;
 737         GLuint uniform_diff_flow_tex, uniform_base_flow_tex;
 738         GLuint uniform_beta_0_tex;
 739         GLuint uniform_smoothness_x_tex, uniform_smoothness_y_tex;
 740         GLuint uniform_gamma, uniform_delta;
 741 };
 742
 743 SetupEquations::SetupEquations()
 744 {
 745         equations_vs_obj = compile_shader(read_file("vs.vert"), GL_VERTEX_SHADER);
 746         equations_fs_obj = compile_shader(read_file("equations.frag"), GL_FRAGMENT_SHADER);
 747         equations_program = link_program(equations_vs_obj, equations_fs_obj);
 748
 749         // Set up the VAO containing all the required position/texcoord data.
 750         glCreateVertexArrays(1, &equations_vao);
 751         glBindVertexArray(equations_vao);
 752         glBindBuffer(GL_ARRAY_BUFFER, vertex_vbo);
 753
 754         GLint position_attrib = glGetAttribLocation(equations_program, "position");
 755         glEnableVertexArrayAttrib(equations_vao, position_attrib);
 756         glVertexAttribPointer(position_attrib, 2, GL_FLOAT, GL_FALSE, 0, BUFFER_OFFSET(0));
 757
 758         uniform_I_x_y_tex = glGetUniformLocation(equations_program, "I_x_y_tex");
 759         uniform_I_t_tex = glGetUniformLocation(equations_program, "I_t_tex");
 760         uniform_diff_flow_tex = glGetUniformLocation(equations_program, "diff_flow_tex");
 761         uniform_base_flow_tex = glGetUniformLocation(equations_program, "base_flow_tex");
 762         uniform_beta_0_tex = glGetUniformLocation(equations_program, "beta_0_tex");
 763         uniform_smoothness_x_tex = glGetUniformLocation(equations_program, "smoothness_x_tex");
 764         uniform_smoothness_y_tex = glGetUniformLocation(equations_program, "smoothness_y_tex");
 765         uniform_gamma = glGetUniformLocation(equations_program, "gamma");
 766         uniform_delta = glGetUniformLocation(equations_program, "delta");
 767 }
 768
 769 void SetupEquations::exec(GLuint I_x_y_tex, GLuint I_t_tex, GLuint diff_flow_tex, GLuint base_flow_tex, GLuint beta_0_tex, GLuint smoothness_x_tex, GLuint smoothness_y_tex, GLuint equation_tex, int level_width, int level_height)
 770 {
 771         glUseProgram(equations_program);
 772
 773         bind_sampler(equations_program, uniform_I_x_y_tex, 0, I_x_y_tex, nearest_sampler);
 774         bind_sampler(equations_program, uniform_I_t_tex, 1, I_t_tex, nearest_sampler);
 775         bind_sampler(equations_program, uniform_diff_flow_tex, 2, diff_flow_tex, nearest_sampler);
 776         bind_sampler(equations_program, uniform_base_flow_tex, 3, base_flow_tex, nearest_sampler);
 777         bind_sampler(equations_program, uniform_beta_0_tex, 4, beta_0_tex, nearest_sampler);
 778         bind_sampler(equations_program, uniform_smoothness_x_tex, 5, smoothness_x_tex, zero_border_sampler);
 779         bind_sampler(equations_program, uniform_smoothness_y_tex, 6, smoothness_y_tex, zero_border_sampler);
 780         glProgramUniform1f(equations_program, uniform_delta, vr_delta);
 781         glProgramUniform1f(equations_program, uniform_gamma, vr_gamma);
 782
 783         glViewport(0, 0, level_width, level_height);
 784         glDisable(GL_BLEND);
 785         glBindVertexArray(equations_vao);
 786         fbos.render_to(equation_tex);
 787         glDrawArrays(GL_TRIANGLE_STRIP, 0, 4);
 788 }
 789
 790 // Actually solve the equation sets made by SetupEquations, by means of
 791 // successive over-relaxation (SOR).
 792 //
 793 // See variational_refinement.txt for more information.
 794 class SOR {
 795 public:
 796         SOR();
 797         void exec(GLuint diff_flow_tex, GLuint equation_tex, GLuint smoothness_x_tex, GLuint smoothness_y_tex, int level_width, int level_height, int num_iterations);
 798
 799 private:
 800         PersistentFBOSet<1> fbos;
 801
 802         GLuint sor_vs_obj;
 803         GLuint sor_fs_obj;
 804         GLuint sor_program;
 805         GLuint sor_vao;
 806
 807         GLuint uniform_diff_flow_tex;
 808         GLuint uniform_equation_tex;
 809         GLuint uniform_smoothness_x_tex, uniform_smoothness_y_tex;
 810         GLuint uniform_phase;
 811 };
 812
 813 SOR::SOR()
 814 {
 815         sor_vs_obj = compile_shader(read_file("sor.vert"), GL_VERTEX_SHADER);
 816         sor_fs_obj = compile_shader(read_file("sor.frag"), GL_FRAGMENT_SHADER);
 817         sor_program = link_program(sor_vs_obj, sor_fs_obj);
 818
 819         // Set up the VAO containing all the required position/texcoord data.
 820         glCreateVertexArrays(1, &sor_vao);
 821         glBindVertexArray(sor_vao);
 822         glBindBuffer(GL_ARRAY_BUFFER, vertex_vbo);
 823
 824         GLint position_attrib = glGetAttribLocation(sor_program, "position");
 825         glEnableVertexArrayAttrib(sor_vao, position_attrib);
 826         glVertexAttribPointer(position_attrib, 2, GL_FLOAT, GL_FALSE, 0, BUFFER_OFFSET(0));
 827
 828         uniform_diff_flow_tex = glGetUniformLocation(sor_program, "diff_flow_tex");
 829         uniform_equation_tex = glGetUniformLocation(sor_program, "equation_tex");
 830         uniform_smoothness_x_tex = glGetUniformLocation(sor_program, "smoothness_x_tex");
 831         uniform_smoothness_y_tex = glGetUniformLocation(sor_program, "smoothness_y_tex");
 832         uniform_phase = glGetUniformLocation(sor_program, "phase");
 833 }
 834
 835 void SOR::exec(GLuint diff_flow_tex, GLuint equation_tex, GLuint smoothness_x_tex, GLuint smoothness_y_tex, int level_width, int level_height, int num_iterations)
 836 {
 837         glUseProgram(sor_program);
 838
 839         bind_sampler(sor_program, uniform_diff_flow_tex, 0, diff_flow_tex, nearest_sampler);
 840         bind_sampler(sor_program, uniform_smoothness_x_tex, 1, smoothness_x_tex, zero_border_sampler);
 841         bind_sampler(sor_program, uniform_smoothness_y_tex, 2, smoothness_y_tex, zero_border_sampler);
 842         bind_sampler(sor_program, uniform_equation_tex, 3, equation_tex, nearest_sampler);
 843
 844         // NOTE: We bind to the texture we are rendering from, but we never write any value
 845         // that we read in the same shader pass (we call discard for red values when we compute
 846         // black, and vice versa), and we have barriers between the passes, so we're fine
 847         // as per the spec.
 848         glViewport(0, 0, level_width, level_height);
 849         glDisable(GL_BLEND);
 850         glBindVertexArray(sor_vao);
 851         fbos.render_to(diff_flow_tex);
 852
 853         for (int i = 0; i < num_iterations; ++i) {
 854                 glProgramUniform1i(sor_program, uniform_phase, 0);
 855                 glDrawArrays(GL_TRIANGLE_STRIP, 0, 4);
 856                 glTextureBarrier();
 857                 glProgramUniform1i(sor_program, uniform_phase, 1);
 858                 glDrawArrays(GL_TRIANGLE_STRIP, 0, 4);
 859                 if (i != num_iterations - 1) {
 860                         glTextureBarrier();
 861                 }
 862         }
 863 }
 864
 865 // Simply add the differential flow found by the variational refinement to the base flow.
 866 // The output is in base_flow_tex; we don't need to make a new texture.
 867 class AddBaseFlow {
 868 public:
 869         AddBaseFlow();
 870         void exec(GLuint base_flow_tex, GLuint diff_flow_tex, int level_width, int level_height);
 871
 872 private:
 873         PersistentFBOSet<1> fbos;
 874
 875         GLuint add_flow_vs_obj;
 876         GLuint add_flow_fs_obj;
 877         GLuint add_flow_program;
 878         GLuint add_flow_vao;
 879
 880         GLuint uniform_diff_flow_tex;
 881 };
 882
 883 AddBaseFlow::AddBaseFlow()
 884 {
 885         add_flow_vs_obj = compile_shader(read_file("vs.vert"), GL_VERTEX_SHADER);
 886         add_flow_fs_obj = compile_shader(read_file("add_base_flow.frag"), GL_FRAGMENT_SHADER);
 887         add_flow_program = link_program(add_flow_vs_obj, add_flow_fs_obj);
 888
 889         // Set up the VAO containing all the required position/texcoord data.
 890         glCreateVertexArrays(1, &add_flow_vao);
 891         glBindVertexArray(add_flow_vao);
 892         glBindBuffer(GL_ARRAY_BUFFER, vertex_vbo);
 893
 894         GLint position_attrib = glGetAttribLocation(add_flow_program, "position");
 895         glEnableVertexArrayAttrib(add_flow_vao, position_attrib);
 896         glVertexAttribPointer(position_attrib, 2, GL_FLOAT, GL_FALSE, 0, BUFFER_OFFSET(0));
 897
 898         uniform_diff_flow_tex = glGetUniformLocation(add_flow_program, "diff_flow_tex");
 899 }
 900
 901 void AddBaseFlow::exec(GLuint base_flow_tex, GLuint diff_flow_tex, int level_width, int level_height)
 902 {
 903         glUseProgram(add_flow_program);
 904
 905         bind_sampler(add_flow_program, uniform_diff_flow_tex, 0, diff_flow_tex, nearest_sampler);
 906
 907         glViewport(0, 0, level_width, level_height);
 908         glEnable(GL_BLEND);
 909         glBlendFunc(GL_ONE, GL_ONE);
 910         glBindVertexArray(add_flow_vao);
 911         fbos.render_to(base_flow_tex);
 912
 913         glDrawArrays(GL_TRIANGLE_STRIP, 0, 4);
 914 }
 915
 916 // Take a copy of the flow, bilinearly interpolated and scaled up.
 917 class ResizeFlow {
 918 public:
 919         ResizeFlow();
 920         void exec(GLuint in_tex, GLuint out_tex, int input_width, int input_height, int output_width, int output_height);
 921
 922 private:
 923         PersistentFBOSet<1> fbos;
 924
 925         GLuint resize_flow_vs_obj;
 926         GLuint resize_flow_fs_obj;
 927         GLuint resize_flow_program;
 928         GLuint resize_flow_vao;
 929
 930         GLuint uniform_flow_tex;
 931         GLuint uniform_scale_factor;
 932 };
 933
 934 ResizeFlow::ResizeFlow()
 935 {
 936         resize_flow_vs_obj = compile_shader(read_file("vs.vert"), GL_VERTEX_SHADER);
 937         resize_flow_fs_obj = compile_shader(read_file("resize_flow.frag"), GL_FRAGMENT_SHADER);
 938         resize_flow_program = link_program(resize_flow_vs_obj, resize_flow_fs_obj);
 939
 940         // Set up the VAO containing all the required position/texcoord data.
 941         glCreateVertexArrays(1, &resize_flow_vao);
 942         glBindVertexArray(resize_flow_vao);
 943         glBindBuffer(GL_ARRAY_BUFFER, vertex_vbo);
 944
 945         GLint position_attrib = glGetAttribLocation(resize_flow_program, "position");
 946         glEnableVertexArrayAttrib(resize_flow_vao, position_attrib);
 947         glVertexAttribPointer(position_attrib, 2, GL_FLOAT, GL_FALSE, 0, BUFFER_OFFSET(0));
 948
 949         uniform_flow_tex = glGetUniformLocation(resize_flow_program, "flow_tex");
 950         uniform_scale_factor = glGetUniformLocation(resize_flow_program, "scale_factor");
 951 }
 952
 953 void ResizeFlow::exec(GLuint flow_tex, GLuint out_tex, int input_width, int input_height, int output_width, int output_height)
 954 {
 955         glUseProgram(resize_flow_program);
 956
 957         bind_sampler(resize_flow_program, uniform_flow_tex, 0, flow_tex, nearest_sampler);
 958
 959         glProgramUniform2f(resize_flow_program, uniform_scale_factor, float(output_width) / input_width, float(output_height) / input_height);
 960
 961         glViewport(0, 0, output_width, output_height);
 962         glDisable(GL_BLEND);
 963         glBindVertexArray(resize_flow_vao);
 964         fbos.render_to(out_tex);
 965
 966         glDrawArrays(GL_TRIANGLE_STRIP, 0, 4);
 967 }
 968
 969 class GPUTimers {
 970 public:
 971         void print();
 972         pair<GLuint, GLuint> begin_timer(const string &name, int level);
 973
 974 private:
 975         struct Timer {
 976                 string name;
 977                 int level;
 978                 pair<GLuint, GLuint> query;
 979         };
 980         vector<Timer> timers;
 981 };
 982
 983 pair<GLuint, GLuint> GPUTimers::begin_timer(const string &name, int level)
 984 {
 985         if (!enable_timing) {
 986                 return make_pair(0, 0);
 987         }
 988
 989         GLuint queries[2];
 990         glGenQueries(2, queries);
 991         glQueryCounter(queries[0], GL_TIMESTAMP);
 992
 993         Timer timer;
 994         timer.name = name;
 995         timer.level = level;
 996         timer.query.first = queries[0];
 997         timer.query.second = queries[1];
 998         timers.push_back(timer);
 999         return timer.query;
1000 }
1001
1002 void GPUTimers::print()
1003 {
1004         for (const Timer &timer : timers) {
1005                 // NOTE: This makes the CPU wait for the GPU.
1006                 GLuint64 time_start, time_end;
1007                 glGetQueryObjectui64v(timer.query.first, GL_QUERY_RESULT, &time_start);
1008                 glGetQueryObjectui64v(timer.query.second, GL_QUERY_RESULT, &time_end);
1009                 //fprintf(stderr, "GPU time used = %.1f ms\n", time_elapsed / 1e6);
1010                 for (int i = 0; i < timer.level * 2; ++i) {
1011                         fprintf(stderr, " ");
1012                 }
1013                 fprintf(stderr, "%-30s %4.1f ms\n", timer.name.c_str(), GLint64(time_end - time_start) / 1e6);
1014         }
1015 }
1016
1017 // A simple RAII class for timing until the end of the scope.
1018 class ScopedTimer {
1019 public:
1020         ScopedTimer(const string &name, GPUTimers *timers)
1021                 : timers(timers), level(0)
1022         {
1023                 query = timers->begin_timer(name, level);
1024         }
1025
1026         ScopedTimer(const string &name, ScopedTimer *parent_timer)
1027                 : timers(parent_timer->timers),
1028                   level(parent_timer->level + 1)
1029         {
1030                 query = timers->begin_timer(name, level);
1031         }
1032
1033         ~ScopedTimer()
1034         {
1035                 end();
1036         }
1037
1038         void end()
1039         {
1040                 if (enable_timing && !ended) {
1041                         glQueryCounter(query.second, GL_TIMESTAMP);
1042                         ended = true;
1043                 }
1044         }
1045
1046 private:
1047         GPUTimers *timers;
1048         int level;
1049         pair<GLuint, GLuint> query;
1050         bool ended = false;
1051 };
1052
1053 class TexturePool {
1054 public:
1055         GLuint get_texture(GLenum format, GLuint width, GLuint height);
1056         void release_texture(GLuint tex_num);
1057
1058 private:
1059         struct Texture {
1060                 GLuint tex_num;
1061                 GLenum format;
1062                 GLuint width, height;
1063                 bool in_use = false;
1064         };
1065         vector<Texture> textures;
1066 };
1067
1068 class DISComputeFlow {
1069 public:
1070         DISComputeFlow(int width, int height);
1071
1072         // Returns a texture that must be released with release_texture()
1073         // after use.
1074         GLuint exec(GLuint tex0, GLuint tex1);
1075
1076         void release_texture(GLuint tex) {
1077                 pool.release_texture(tex);
1078         }
1079
1080 private:
1081         int width, height;
1082         GLuint initial_flow_tex;
1083         TexturePool pool;
1084
1085         // The various passes.
1086         Sobel sobel;
1087         MotionSearch motion_search;
1088         Densify densify;
1089         Prewarp prewarp;
1090         Derivatives derivatives;
1091         ComputeSmoothness compute_smoothness;
1092         SetupEquations setup_equations;
1093         SOR sor;
1094         AddBaseFlow add_base_flow;
1095         ResizeFlow resize_flow;
1096 };
1097
1098 DISComputeFlow::DISComputeFlow(int width, int height)
1099         : width(width), height(height)
1100 {
1101         // Make some samplers.
1102         glCreateSamplers(1, &nearest_sampler);
1103         glSamplerParameteri(nearest_sampler, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
1104         glSamplerParameteri(nearest_sampler, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
1105         glSamplerParameteri(nearest_sampler, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
1106         glSamplerParameteri(nearest_sampler, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
1107
1108         glCreateSamplers(1, &linear_sampler);
1109         glSamplerParameteri(linear_sampler, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
1110         glSamplerParameteri(linear_sampler, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
1111         glSamplerParameteri(linear_sampler, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
1112         glSamplerParameteri(linear_sampler, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
1113
1114         // The smoothness is sampled so that once we get to a smoothness involving
1115         // a value outside the border, the diffusivity between the two becomes zero.
1116         // Similarly, gradients are zero outside the border, since the edge is taken
1117         // to be constant.
1118         glCreateSamplers(1, &zero_border_sampler);
1119         glSamplerParameteri(zero_border_sampler, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
1120         glSamplerParameteri(zero_border_sampler, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
1121         glSamplerParameteri(zero_border_sampler, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_BORDER);
1122         glSamplerParameteri(zero_border_sampler, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_BORDER);
1123         float zero[] = { 0.0f, 0.0f, 0.0f, 0.0f };
1124         glSamplerParameterfv(zero_border_sampler, GL_TEXTURE_BORDER_COLOR, zero);
1125
1126         // Initial flow is zero, 1x1.
1127         glCreateTextures(GL_TEXTURE_2D, 1, &initial_flow_tex);
1128         glTextureStorage2D(initial_flow_tex, 1, GL_RG16F, 1, 1);
1129         glClearTexImage(initial_flow_tex, 0, GL_RG, GL_FLOAT, nullptr);
1130 }
1131
1132 GLuint DISComputeFlow::exec(GLuint tex0, GLuint tex1)
1133 {
1134         int prev_level_width = 1, prev_level_height = 1;
1135         GLuint prev_level_flow_tex = initial_flow_tex;
1136
1137         GPUTimers timers;
1138
1139         ScopedTimer total_timer("Total", &timers);
1140         for (int level = coarsest_level; level >= int(finest_level); --level) {
1141                 char timer_name[256];
1142                 snprintf(timer_name, sizeof(timer_name), "Level %d", level);
1143                 ScopedTimer level_timer(timer_name, &total_timer);
1144
1145                 int level_width = width >> level;
1146                 int level_height = height >> level;
1147                 float patch_spacing_pixels = patch_size_pixels * (1.0f - patch_overlap_ratio);
1148
1149                 // Make sure we have patches at least every Nth pixel, e.g. for width=9
1150                 // and patch_spacing=3 (the default), we put out patch centers in
1151                 // x=0, x=3, x=6, x=9, which is four patches. The fragment shader will
1152                 // lock all the centers to integer coordinates if needed.
1153                 int width_patches = 1 + ceil(level_width / patch_spacing_pixels);
1154                 int height_patches = 1 + ceil(level_height / patch_spacing_pixels);
1155
1156                 // Make sure we always read from the correct level; the chosen
1157                 // mipmapping could otherwise be rather unpredictable, especially
1158                 // during motion search.
1159                 GLuint tex0_view, tex1_view;
1160                 glGenTextures(1, &tex0_view);
1161                 glTextureView(tex0_view, GL_TEXTURE_2D, tex0, GL_R8, level, 1, 0, 1);
1162                 glGenTextures(1, &tex1_view);
1163                 glTextureView(tex1_view, GL_TEXTURE_2D, tex1, GL_R8, level, 1, 0, 1);
1164
1165                 // Create a new texture; we could be fancy and render use a multi-level
1166                 // texture, but meh.
1167                 GLuint grad0_tex = pool.get_texture(GL_RG16F, level_width, level_height);
1168
1169                 // Find the derivative.
1170                 {
1171                         ScopedTimer timer("Sobel", &level_timer);
1172                         sobel.exec(tex0_view, grad0_tex, level_width, level_height);
1173                 }
1174
1175                 // Motion search to find the initial flow. We use the flow from the previous
1176                 // level (sampled bilinearly; no fancy tricks) as a guide, then search from there.
1177
1178                 // Create an output flow texture.
1179                 GLuint flow_out_tex = pool.get_texture(GL_RGB16F, width_patches, height_patches);
1180
1181                 // And draw.
1182                 {
1183                         ScopedTimer timer("Motion search", &level_timer);
1184                         motion_search.exec(tex0_view, tex1_view, grad0_tex, prev_level_flow_tex, flow_out_tex, level_width, level_height, prev_level_width, prev_level_height, width_patches, height_patches);
1185                 }
1186                 pool.release_texture(grad0_tex);
1187
1188                 // Densification.
1189
1190                 // Set up an output texture (initially zero).
1191                 GLuint dense_flow_tex = pool.get_texture(GL_RGB16F, level_width, level_height);
1192                 glClearTexImage(dense_flow_tex, 0, GL_RGB, GL_FLOAT, nullptr);
1193
1194                 // And draw.
1195                 {
1196                         ScopedTimer timer("Densification", &level_timer);
1197                         densify.exec(tex0_view, tex1_view, flow_out_tex, dense_flow_tex, level_width, level_height, width_patches, height_patches);
1198                 }
1199                 pool.release_texture(flow_out_tex);
1200
1201                 // Everything below here in the loop belongs to variational refinement.
1202                 ScopedTimer varref_timer("Variational refinement", &level_timer);
1203
1204                 // Prewarping; create I and I_t, and a normalized base flow (so we don't
1205                 // have to normalize it over and over again, and also save some bandwidth).
1206                 //
1207                 // During the entire rest of the variational refinement, flow will be measured
1208                 // in pixels, not 0..1 normalized OpenGL texture coordinates.
1209                 // This is because variational refinement depends so heavily on derivatives,
1210                 // which are measured in intensity levels per pixel.
1211                 GLuint I_tex = pool.get_texture(GL_R16F, level_width, level_height);
1212                 GLuint I_t_tex = pool.get_texture(GL_R16F, level_width, level_height);
1213                 GLuint base_flow_tex = pool.get_texture(GL_RG16F, level_width, level_height);
1214                 {
1215                         ScopedTimer timer("Prewarping", &varref_timer);
1216                         prewarp.exec(tex0_view, tex1_view, dense_flow_tex, I_tex, I_t_tex, base_flow_tex, level_width, level_height);
1217                 }
1218                 pool.release_texture(dense_flow_tex);
1219                 glDeleteTextures(1, &tex0_view);
1220                 glDeleteTextures(1, &tex1_view);
1221
1222                 // Calculate I_x and I_y. We're only calculating first derivatives;
1223                 // the others will be taken on-the-fly in order to sample from fewer
1224                 // textures overall, since sampling from the L1 cache is cheap.
1225                 // (TODO: Verify that this is indeed faster than making separate
1226                 // double-derivative textures.)
1227                 GLuint I_x_y_tex = pool.get_texture(GL_RG16F, level_width, level_height);
1228                 GLuint beta_0_tex = pool.get_texture(GL_R16F, level_width, level_height);
1229                 {
1230                         ScopedTimer timer("First derivatives", &varref_timer);
1231                         derivatives.exec(I_tex, I_x_y_tex, beta_0_tex, level_width, level_height);
1232                 }
1233                 pool.release_texture(I_tex);
1234
1235                 // We need somewhere to store du and dv (the flow increment, relative
1236                 // to the non-refined base flow u0 and v0). It starts at zero.
1237                 GLuint du_dv_tex = pool.get_texture(GL_RG16F, level_width, level_height);
1238                 glClearTexImage(du_dv_tex, 0, GL_RG, GL_FLOAT, nullptr);
1239
1240                 // And for smoothness.
1241                 GLuint smoothness_x_tex = pool.get_texture(GL_R16F, level_width, level_height);
1242                 GLuint smoothness_y_tex = pool.get_texture(GL_R16F, level_width, level_height);
1243
1244                 // And finally for the equation set. See SetupEquations for
1245                 // the storage format.
1246                 GLuint equation_tex = pool.get_texture(GL_RGBA32UI, level_width, level_height);
1247
1248                 for (int outer_idx = 0; outer_idx < level + 1; ++outer_idx) {
1249                         // Calculate the smoothness terms between the neighboring pixels,
1250                         // both in x and y direction.
1251                         {
1252                                 ScopedTimer timer("Compute smoothness", &varref_timer);
1253                                 compute_smoothness.exec(base_flow_tex, du_dv_tex, smoothness_x_tex, smoothness_y_tex, level_width, level_height);
1254                         }
1255
1256                         // Set up the 2x2 equation system for each pixel.
1257                         {
1258                                 ScopedTimer timer("Set up equations", &varref_timer);
1259                                 setup_equations.exec(I_x_y_tex, I_t_tex, du_dv_tex, base_flow_tex, beta_0_tex, smoothness_x_tex, smoothness_y_tex, equation_tex, level_width, level_height);
1260                         }
1261
1262                         // Run a few SOR (or quasi-SOR, since we're not really Jacobi) iterations.
1263                         // Note that these are to/from the same texture.
1264                         {
1265                                 ScopedTimer timer("SOR", &varref_timer);
1266                                 sor.exec(du_dv_tex, equation_tex, smoothness_x_tex, smoothness_y_tex, level_width, level_height, 5);
1267                         }
1268                 }
1269
1270                 pool.release_texture(I_t_tex);
1271                 pool.release_texture(I_x_y_tex);
1272                 pool.release_texture(beta_0_tex);
1273                 pool.release_texture(smoothness_x_tex);
1274                 pool.release_texture(smoothness_y_tex);
1275                 pool.release_texture(equation_tex);
1276
1277                 // Add the differential flow found by the variational refinement to the base flow,
1278                 // giving the final flow estimate for this level.
1279                 // The output is in diff_flow_tex; we don't need to make a new texture.
1280                 //
1281                 // Disabling this doesn't save any time (although we could easily make it so that
1282                 // it is more efficient), but it helps debug the motion search.
1283                 if (enable_variational_refinement) {
1284                         ScopedTimer timer("Add differential flow", &varref_timer);
1285                         add_base_flow.exec(base_flow_tex, du_dv_tex, level_width, level_height);
1286                 }
1287                 pool.release_texture(du_dv_tex);
1288
1289                 if (prev_level_flow_tex != initial_flow_tex) {
1290                         pool.release_texture(prev_level_flow_tex);
1291                 }
1292                 prev_level_flow_tex = base_flow_tex;
1293                 prev_level_width = level_width;
1294                 prev_level_height = level_height;
1295         }
1296         total_timer.end();
1297
1298         timers.print();
1299
1300         // Scale up the flow to the final size (if needed).
1301         if (finest_level == 0) {
1302                 return prev_level_flow_tex;
1303         } else {
1304                 GLuint final_tex = pool.get_texture(GL_RG16F, width, height);
1305                 resize_flow.exec(prev_level_flow_tex, final_tex, prev_level_width, prev_level_height, width, height);
1306                 pool.release_texture(prev_level_flow_tex);
1307                 return final_tex;
1308         }
1309 }
1310
1311 GLuint TexturePool::get_texture(GLenum format, GLuint width, GLuint height)
1312 {
1313         for (Texture &tex : textures) {
1314                 if (!tex.in_use && tex.format == format &&
1315                     tex.width == width && tex.height == height) {
1316                         tex.in_use = true;
1317                         return tex.tex_num;
1318                 }
1319         }
1320
1321         Texture tex;
1322         glCreateTextures(GL_TEXTURE_2D, 1, &tex.tex_num);
1323         glTextureStorage2D(tex.tex_num, 1, format, width, height);
1324         tex.format = format;
1325         tex.width = width;
1326         tex.height = height;
1327         tex.in_use = true;
1328         textures.push_back(tex);
1329         return tex.tex_num;
1330 }
1331
1332 void TexturePool::release_texture(GLuint tex_num)
1333 {
1334         for (Texture &tex : textures) {
1335                 if (tex.tex_num == tex_num) {
1336                         assert(tex.in_use);
1337                         tex.in_use = false;
1338                         return;
1339                 }
1340         }
1341         assert(false);
1342 }
1343
1344 // OpenGL uses a bottom-left coordinate system, .flo files use a top-left coordinate system.
1345 void flip_coordinate_system(float *dense_flow, unsigned width, unsigned height)
1346 {
1347         for (unsigned i = 0; i < width * height; ++i) {
1348                 dense_flow[i * 2 + 1] = -dense_flow[i * 2 + 1];
1349         }
1350 }
1351
1352 void write_flow(const char *filename, const float *dense_flow, unsigned width, unsigned height)
1353 {
1354         FILE *flowfp = fopen(filename, "wb");
1355         fprintf(flowfp, "FEIH");
1356         fwrite(&width, 4, 1, flowfp);
1357         fwrite(&height, 4, 1, flowfp);
1358         for (unsigned y = 0; y < height; ++y) {
1359                 int yy = height - y - 1;
1360                 fwrite(&dense_flow[yy * width * 2], width * 2 * sizeof(float), 1, flowfp);
1361         }
1362         fclose(flowfp);
1363 }
1364
1365 void write_ppm(const char *filename, const float *dense_flow, unsigned width, unsigned height)
1366 {
1367         FILE *fp = fopen(filename, "wb");
1368         fprintf(fp, "P6\n%d %d\n255\n", width, height);
1369         for (unsigned y = 0; y < unsigned(height); ++y) {
1370                 int yy = height - y - 1;
1371                 for (unsigned x = 0; x < unsigned(width); ++x) {
1372                         float du = dense_flow[(yy * width + x) * 2 + 0];
1373                         float dv = dense_flow[(yy * width + x) * 2 + 1];
1374
1375                         uint8_t r, g, b;
1376                         flow2rgb(du, dv, &r, &g, &b);
1377                         putc(r, fp);
1378                         putc(g, fp);
1379                         putc(b, fp);
1380                 }
1381         }
1382         fclose(fp);
1383 }
1384
1385 void finish_one_read(GLuint width, GLuint height)
1386 {
1387         assert(!reads_in_progress.empty());
1388         ReadInProgress read = reads_in_progress.front();
1389         reads_in_progress.pop_front();
1390
1391         unique_ptr<float[]> flow(new float[width * height * 2]);
1392         void *buf = glMapNamedBufferRange(read.pbo, 0, width * height * 2 * sizeof(float), GL_MAP_READ_BIT);  // Blocks if the read isn't done yet.
1393         memcpy(flow.get(), buf, width * height * 2 * sizeof(float));
1394         glUnmapNamedBuffer(read.pbo);
1395         spare_pbos.push(read.pbo);
1396
1397         flip_coordinate_system(flow.get(), width, height);
1398         if (!read.flow_filename.empty()) {
1399                 write_flow(read.flow_filename.c_str(), flow.get(), width, height);
1400                 fprintf(stderr, "%s %s -> %s\n", read.filename0.c_str(), read.filename1.c_str(), read.flow_filename.c_str());
1401         }
1402         if (!read.ppm_filename.empty()) {
1403                 write_ppm(read.ppm_filename.c_str(), flow.get(), width, height);
1404         }
1405 }
1406
1407 void schedule_read(GLuint tex, GLuint width, GLuint height, const char *filename0, const char *filename1, const char *flow_filename, const char *ppm_filename)
1408 {
1409         if (spare_pbos.empty()) {
1410                 finish_one_read(width, height);
1411         }
1412         assert(!spare_pbos.empty());
1413         reads_in_progress.emplace_back(ReadInProgress{ spare_pbos.top(), filename0, filename1, flow_filename, ppm_filename });
1414         glBindBuffer(GL_PIXEL_PACK_BUFFER, spare_pbos.top());
1415         spare_pbos.pop();
1416         glGetTextureImage(tex, 0, GL_RG, GL_FLOAT, width * height * 2 * sizeof(float), nullptr);
1417         glBindBuffer(GL_PIXEL_PACK_BUFFER, 0);
1418 }
1419
1420 int main(int argc, char **argv)
1421 {
1422         static const option long_options[] = {
1423                 { "smoothness-relative-weight", required_argument, 0, 's' },  // alpha.
1424                 { "intensity-relative-weight", required_argument, 0, 'i' },  // delta.
1425                 { "gradient-relative-weight", required_argument, 0, 'g' },  // gamma.
1426                 { "disable-timing", no_argument, 0, 1000 },
1427                 { "ignore-variational-refinement", no_argument, 0, 1001 }  // Still calculates it, just doesn't apply it.
1428         };
1429
1430         for ( ;; ) {
1431                 int option_index = 0;
1432                 int c = getopt_long(argc, argv, "s:i:g:", long_options, &option_index);
1433
1434                 if (c == -1) {
1435                         break;
1436                 }
1437                 switch (c) {
1438                 case 's':
1439                         vr_alpha = atof(optarg);
1440                         break;
1441                 case 'i':
1442                         vr_delta = atof(optarg);
1443                         break;
1444                 case 'g':
1445                         vr_gamma = atof(optarg);
1446                         break;
1447                 case 1000:
1448                         enable_timing = false;
1449                         break;
1450                 case 1001:
1451                         enable_variational_refinement = false;
1452                         break;
1453                 default:
1454                         fprintf(stderr, "Unknown option '%s'\n", argv[option_index]);
1455                         exit(1);
1456                 };
1457         }
1458
1459         if (SDL_Init(SDL_INIT_EVERYTHING) == -1) {
1460                 fprintf(stderr, "SDL_Init failed: %s\n", SDL_GetError());
1461                 exit(1);
1462         }
1463         SDL_GL_SetAttribute(SDL_GL_ALPHA_SIZE, 8);
1464         SDL_GL_SetAttribute(SDL_GL_DEPTH_SIZE, 0);
1465         SDL_GL_SetAttribute(SDL_GL_STENCIL_SIZE, 0);
1466         SDL_GL_SetAttribute(SDL_GL_DOUBLEBUFFER, 1);
1467
1468         SDL_GL_SetAttribute(SDL_GL_CONTEXT_PROFILE_MASK, SDL_GL_CONTEXT_PROFILE_CORE);
1469         SDL_GL_SetAttribute(SDL_GL_CONTEXT_MAJOR_VERSION, 4);
1470         SDL_GL_SetAttribute(SDL_GL_CONTEXT_MINOR_VERSION, 5);
1471         // SDL_GL_SetAttribute(SDL_GL_CONTEXT_FLAGS, SDL_GL_CONTEXT_DEBUG_FLAG);
1472         SDL_Window *window = SDL_CreateWindow("OpenGL window",
1473                 SDL_WINDOWPOS_UNDEFINED,
1474                 SDL_WINDOWPOS_UNDEFINED,
1475                 64, 64,
1476                 SDL_WINDOW_OPENGL | SDL_WINDOW_HIDDEN);
1477         SDL_GLContext context = SDL_GL_CreateContext(window);
1478         assert(context != nullptr);
1479
1480         const char *filename0 = argc >= (optind + 1) ? argv[optind] : "test1499.png";
1481         const char *filename1 = argc >= (optind + 2) ? argv[optind + 1] : "test1500.png";
1482         const char *flow_filename = argc >= (optind + 3) ? argv[optind + 2] : "flow.flo";
1483
1484         // Load pictures.
1485         unsigned width1, height1, width2, height2;
1486         GLuint tex0 = load_texture(filename0, &width1, &height1);
1487         GLuint tex1 = load_texture(filename1, &width2, &height2);
1488
1489         if (width1 != width2 || height1 != height2) {
1490                 fprintf(stderr, "Image dimensions don't match (%dx%d versus %dx%d)\n",
1491                         width1, height1, width2, height2);
1492                 exit(1);
1493         }
1494
1495         // Set up some PBOs to do asynchronous readback.
1496         GLuint pbos[5];
1497         glCreateBuffers(5, pbos);
1498         for (int i = 0; i < 5; ++i) {
1499                 glNamedBufferData(pbos[i], width1 * height1 * 2 * sizeof(float), nullptr, GL_STREAM_READ);
1500                 spare_pbos.push(pbos[i]);
1501         }
1502
1503         // FIXME: Should be part of DISComputeFlow (but needs to be initialized
1504         // before all the render passes).
1505         float vertices[] = {
1506                 0.0f, 1.0f,
1507                 0.0f, 0.0f,
1508                 1.0f, 1.0f,
1509                 1.0f, 0.0f,
1510         };
1511         glCreateBuffers(1, &vertex_vbo);
1512         glNamedBufferData(vertex_vbo, sizeof(vertices), vertices, GL_STATIC_DRAW);
1513         glBindBuffer(GL_ARRAY_BUFFER, vertex_vbo);
1514
1515         int levels = find_num_levels(width1, height1);
1516         GLuint tex0_gray, tex1_gray;
1517         glCreateTextures(GL_TEXTURE_2D, 1, &tex0_gray);
1518         glCreateTextures(GL_TEXTURE_2D, 1, &tex1_gray);
1519         glTextureStorage2D(tex0_gray, levels, GL_R8, width1, height1);
1520         glTextureStorage2D(tex1_gray, levels, GL_R8, width1, height1);
1521
1522         GrayscaleConversion gray;
1523         gray.exec(tex0, tex0_gray, width1, height1);
1524         glDeleteTextures(1, &tex0);
1525         glGenerateTextureMipmap(tex0_gray);
1526
1527         gray.exec(tex1, tex1_gray, width1, height1);
1528         glDeleteTextures(1, &tex1);
1529         glGenerateTextureMipmap(tex1_gray);
1530
1531         DISComputeFlow compute_flow(width1, height1);
1532         GLuint final_tex = compute_flow.exec(tex0_gray, tex1_gray);
1533
1534         schedule_read(final_tex, width1, height1, filename0, filename1, flow_filename, "flow.ppm");
1535         compute_flow.release_texture(final_tex);
1536
1537         // See if there are more flows on the command line (ie., more than three arguments),
1538         // and if so, process them.
1539         int num_flows = (argc - optind) / 3;
1540         for (int i = 1; i < num_flows; ++i) {
1541                 const char *filename0 = argv[optind + i * 3 + 0];
1542                 const char *filename1 = argv[optind + i * 3 + 1];
1543                 const char *flow_filename = argv[optind + i * 3 + 2];
1544                 GLuint width, height;
1545                 GLuint tex0 = load_texture(filename0, &width, &height);
1546                 if (width != width1 || height != height1) {
1547                         fprintf(stderr, "%s: Image dimensions don't match (%dx%d versus %dx%d)\n",
1548                                 filename0, width, height, width1, height1);
1549                         exit(1);
1550                 }
1551                 gray.exec(tex0, tex0_gray, width, height);
1552                 glGenerateTextureMipmap(tex0_gray);
1553                 glDeleteTextures(1, &tex0);
1554
1555                 GLuint tex1 = load_texture(filename1, &width, &height);
1556                 if (width != width1 || height != height1) {
1557                         fprintf(stderr, "%s: Image dimensions don't match (%dx%d versus %dx%d)\n",
1558                                 filename1, width, height, width1, height1);
1559                         exit(1);
1560                 }
1561                 gray.exec(tex1, tex1_gray, width, height);
1562                 glGenerateTextureMipmap(tex1_gray);
1563                 glDeleteTextures(1, &tex1);
1564
1565                 GLuint final_tex = compute_flow.exec(tex0_gray, tex1_gray);
1566
1567                 schedule_read(final_tex, width1, height1, filename0, filename1, flow_filename, "");
1568                 compute_flow.release_texture(final_tex);
1569         }
1570         glDeleteTextures(1, &tex0_gray);
1571         glDeleteTextures(1, &tex1_gray);
1572
1573         while (!reads_in_progress.empty()) {
1574                 finish_one_read(width1, height1);
1575         }
1576 }