]> git.sesse.net Git - nageru/blob - flow.cpp
dba6607b0eb3e967fcb2786bd31c5051ad72797a
[nageru] / flow.cpp
1 #define NO_SDL_GLEXT 1
2
3 #define WIDTH 1280
4 #define HEIGHT 720
5
6 #include <epoxy/gl.h>
7
8 #include <SDL2/SDL.h>
9 #include <SDL2/SDL_error.h>
10 #include <SDL2/SDL_events.h>
11 #include <SDL2/SDL_image.h>
12 #include <SDL2/SDL_keyboard.h>
13 #include <SDL2/SDL_mouse.h>
14 #include <SDL2/SDL_video.h>
15
16 #include <assert.h>
17 #include <stdio.h>
18
19 #include <algorithm>
20 #include <memory>
21
22 #define BUFFER_OFFSET(i) ((char *)nullptr + (i))
23
24 using namespace std;
25
26 // Operating point 3 (10 Hz on CPU, excluding preprocessing).
27 constexpr float patch_overlap_ratio = 0.75f;
28 constexpr unsigned coarsest_level = 5;
29 constexpr unsigned finest_level = 1;
30 constexpr unsigned patch_size_pixels = 12;
31
32 // Some global OpenGL objects.
33 GLuint nearest_sampler, linear_sampler, mipmap_sampler;
34 GLuint vertex_vbo;
35
36 string read_file(const string &filename)
37 {
38         FILE *fp = fopen(filename.c_str(), "r");
39         if (fp == nullptr) {
40                 perror(filename.c_str());
41                 exit(1);
42         }
43
44         int ret = fseek(fp, 0, SEEK_END);
45         if (ret == -1) {
46                 perror("fseek(SEEK_END)");
47                 exit(1);
48         }
49
50         int size = ftell(fp);
51
52         ret = fseek(fp, 0, SEEK_SET);
53         if (ret == -1) {
54                 perror("fseek(SEEK_SET)");
55                 exit(1);
56         }
57
58         string str;
59         str.resize(size);
60         ret = fread(&str[0], size, 1, fp);
61         if (ret == -1) {
62                 perror("fread");
63                 exit(1);
64         }
65         if (ret == 0) {
66                 fprintf(stderr, "Short read when trying to read %d bytes from %s\n",
67                                 size, filename.c_str());
68                 exit(1);
69         }
70         fclose(fp);
71
72         return str;
73 }
74
75
76 GLuint compile_shader(const string &shader_src, GLenum type)
77 {
78         GLuint obj = glCreateShader(type);
79         const GLchar* source[] = { shader_src.data() };
80         const GLint length[] = { (GLint)shader_src.size() };
81         glShaderSource(obj, 1, source, length);
82         glCompileShader(obj);
83
84         GLchar info_log[4096];
85         GLsizei log_length = sizeof(info_log) - 1;
86         glGetShaderInfoLog(obj, log_length, &log_length, info_log);
87         info_log[log_length] = 0;
88         if (strlen(info_log) > 0) {
89                 fprintf(stderr, "Shader compile log: %s\n", info_log);
90         }
91
92         GLint status;
93         glGetShaderiv(obj, GL_COMPILE_STATUS, &status);
94         if (status == GL_FALSE) {
95                 // Add some line numbers to easier identify compile errors.
96                 string src_with_lines = "/*   1 */ ";
97                 size_t lineno = 1;
98                 for (char ch : shader_src) {
99                         src_with_lines.push_back(ch);
100                         if (ch == '\n') {
101                                 char buf[32];
102                                 snprintf(buf, sizeof(buf), "/* %3zu */ ", ++lineno);
103                                 src_with_lines += buf;
104                         }
105                 }
106
107                 fprintf(stderr, "Failed to compile shader:\n%s\n", src_with_lines.c_str());
108                 exit(1);
109         }
110
111         return obj;
112 }
113
114
115 GLuint load_texture(const char *filename, unsigned width, unsigned height)
116 {
117         FILE *fp = fopen(filename, "rb");
118         if (fp == nullptr) {
119                 perror(filename);
120                 exit(1);
121         }
122         unique_ptr<uint8_t[]> pix(new uint8_t[width * height]);
123         if (fread(pix.get(), width * height, 1, fp) != 1) {
124                 fprintf(stderr, "Short read from %s\n", filename);
125                 exit(1);
126         }
127         fclose(fp);
128
129         // Convert to bottom-left origin.
130         for (unsigned y = 0; y < height / 2; ++y) {
131                 unsigned y2 = height - 1 - y;
132                 swap_ranges(&pix[y * width], &pix[y * width + width], &pix[y2 * width]);
133         }
134
135         int levels = 1;
136         for (int w = width, h = height; w > 1 || h > 1; ) {
137                 w >>= 1;
138                 h >>= 1;
139                 ++levels;
140         }
141
142         GLuint tex;
143         glCreateTextures(GL_TEXTURE_2D, 1, &tex);
144         glTextureStorage2D(tex, levels, GL_R8, width, height);
145         glTextureSubImage2D(tex, 0, 0, 0, width, height, GL_RED, GL_UNSIGNED_BYTE, pix.get());
146         glGenerateTextureMipmap(tex);
147
148         return tex;
149 }
150
151 GLuint link_program(GLuint vs_obj, GLuint fs_obj)
152 {
153         GLuint program = glCreateProgram();
154         glAttachShader(program, vs_obj);
155         glAttachShader(program, fs_obj);
156         glLinkProgram(program);
157         GLint success;
158         glGetProgramiv(program, GL_LINK_STATUS, &success);
159         if (success == GL_FALSE) {
160                 GLchar error_log[1024] = {0};
161                 glGetProgramInfoLog(program, 1024, nullptr, error_log);
162                 fprintf(stderr, "Error linking program: %s\n", error_log);
163                 exit(1);
164         }
165         return program;
166 }
167
168 GLuint generate_vbo(GLint size, GLsizeiptr data_size, const GLvoid *data)
169 {
170         GLuint vbo;
171         glCreateBuffers(1, &vbo);
172         glBufferData(GL_ARRAY_BUFFER, data_size, data, GL_STATIC_DRAW);
173         glNamedBufferData(vbo, data_size, data, GL_STATIC_DRAW);
174         return vbo;
175 }
176
177 GLuint fill_vertex_attribute(GLuint vao, GLuint glsl_program_num, const string &attribute_name, GLint size, GLenum type, GLsizeiptr data_size, const GLvoid *data)
178 {
179         int attrib = glGetAttribLocation(glsl_program_num, attribute_name.c_str());
180         if (attrib == -1) {
181                 return -1;
182         }
183
184         GLuint vbo = generate_vbo(size, data_size, data);
185
186         glBindBuffer(GL_ARRAY_BUFFER, vbo);
187         glEnableVertexArrayAttrib(vao, attrib);
188         glVertexAttribPointer(attrib, size, type, GL_FALSE, 0, BUFFER_OFFSET(0));
189         glBindBuffer(GL_ARRAY_BUFFER, 0);
190
191         return vbo;
192 }
193
194 void bind_sampler(GLuint program, const char *uniform_name, GLuint texture_unit, GLuint tex, GLuint sampler)
195 {
196         GLint location = glGetUniformLocation(program, uniform_name);
197         if (location == -1) {
198                 return;
199         }
200
201         glBindTextureUnit(texture_unit, tex);
202         glBindSampler(texture_unit, sampler);
203         glProgramUniform1i(program, location, texture_unit);
204 }
205
206 // Compute gradients in every point, used for the motion search.
207 // The DIS paper doesn't actually mention how these are computed,
208 // but seemingly, a 3x3 Sobel operator is used here (at least in
209 // later versions of the code), while a [1 -8 0 8 -1] kernel is
210 // used for all the derivatives in the variational refinement part
211 // (which borrows code from DeepFlow). This is inconsistent,
212 // but I guess we're better off with staying with the original
213 // decisions until we actually know having different ones would be better.
214 class Sobel {
215 public:
216         Sobel();
217         void exec(GLint tex0_view, GLint grad0_tex, int level_width, int level_height);
218
219 private:
220         GLuint sobel_vs_obj;
221         GLuint sobel_fs_obj;
222         GLuint sobel_program;
223         GLuint sobel_vao;
224 };
225
226 Sobel::Sobel()
227 {
228         sobel_vs_obj = compile_shader(read_file("vs.vert"), GL_VERTEX_SHADER);
229         sobel_fs_obj = compile_shader(read_file("sobel.frag"), GL_FRAGMENT_SHADER);
230         sobel_program = link_program(sobel_vs_obj, sobel_fs_obj);
231
232         // Set up the VAO containing all the required position/texcoord data.
233         glCreateVertexArrays(1, &sobel_vao);
234         glBindVertexArray(sobel_vao);
235
236         GLint position_attrib = glGetAttribLocation(sobel_program, "position");
237         glEnableVertexArrayAttrib(sobel_vao, position_attrib);
238         glVertexAttribPointer(position_attrib, 2, GL_FLOAT, GL_FALSE, 0, BUFFER_OFFSET(0));
239
240         GLint texcoord_attrib = glGetAttribLocation(sobel_program, "texcoord");
241         glEnableVertexArrayAttrib(sobel_vao, texcoord_attrib);
242         glVertexAttribPointer(texcoord_attrib, 2, GL_FLOAT, GL_FALSE, 0, BUFFER_OFFSET(0));
243 }
244
245 void Sobel::exec(GLint tex0_view, GLint grad0_tex, int level_width, int level_height)
246 {
247         glUseProgram(sobel_program);
248         glBindTextureUnit(0, tex0_view);
249         glBindSampler(0, nearest_sampler);
250         glProgramUniform1i(sobel_program, glGetUniformLocation(sobel_program, "tex"), 0);
251         glProgramUniform1f(sobel_program, glGetUniformLocation(sobel_program, "inv_width"), 1.0f / level_width);
252         glProgramUniform1f(sobel_program, glGetUniformLocation(sobel_program, "inv_height"), 1.0f / level_height);
253
254         GLuint grad0_fbo;  // TODO: cleanup
255         glCreateFramebuffers(1, &grad0_fbo);
256         glNamedFramebufferTexture(grad0_fbo, GL_COLOR_ATTACHMENT0, grad0_tex, 0);
257
258         glViewport(0, 0, level_width, level_height);
259         glBindFramebuffer(GL_FRAMEBUFFER, grad0_fbo);
260         glBindVertexArray(sobel_vao);
261         glUseProgram(sobel_program);
262         glDisable(GL_BLEND);
263         glDrawArrays(GL_TRIANGLE_STRIP, 0, 4);
264 }
265
266 // Motion search to find the initial flow. See motion_search.frag for documentation.
267 class MotionSearch {
268 public:
269         MotionSearch();
270         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 width_patches, int height_patches);
271
272 private:
273         GLuint motion_vs_obj;
274         GLuint motion_fs_obj;
275         GLuint motion_search_program;
276         GLuint motion_search_vao;
277 };
278
279 MotionSearch::MotionSearch()
280 {
281         motion_vs_obj = compile_shader(read_file("motion_search.vert"), GL_VERTEX_SHADER);
282         motion_fs_obj = compile_shader(read_file("motion_search.frag"), GL_FRAGMENT_SHADER);
283         motion_search_program = link_program(motion_vs_obj, motion_fs_obj);
284
285         // Set up the VAO containing all the required position/texcoord data.
286         glCreateVertexArrays(1, &motion_search_vao);
287         glBindVertexArray(motion_search_vao);
288         glBindBuffer(GL_ARRAY_BUFFER, vertex_vbo);
289
290         GLint position_attrib = glGetAttribLocation(motion_search_program, "position");
291         glEnableVertexArrayAttrib(motion_search_vao, position_attrib);
292         glVertexAttribPointer(position_attrib, 2, GL_FLOAT, GL_FALSE, 0, BUFFER_OFFSET(0));
293
294         GLint texcoord_attrib = glGetAttribLocation(motion_search_program, "texcoord");
295         glEnableVertexArrayAttrib(motion_search_vao, texcoord_attrib);
296         glVertexAttribPointer(texcoord_attrib, 2, GL_FLOAT, GL_FALSE, 0, BUFFER_OFFSET(0));
297 }
298
299 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 width_patches, int height_patches)
300 {
301         glUseProgram(motion_search_program);
302
303         bind_sampler(motion_search_program, "image0_tex", 0, tex0_view, nearest_sampler);
304         bind_sampler(motion_search_program, "image1_tex", 1, tex1_view, linear_sampler);
305         bind_sampler(motion_search_program, "grad0_tex", 2, grad0_tex, nearest_sampler);
306         bind_sampler(motion_search_program, "flow_tex", 3, flow_tex, linear_sampler);
307
308         glProgramUniform1f(motion_search_program, glGetUniformLocation(motion_search_program, "image_width"), level_width);
309         glProgramUniform1f(motion_search_program, glGetUniformLocation(motion_search_program, "image_height"), level_height);
310         glProgramUniform1f(motion_search_program, glGetUniformLocation(motion_search_program, "inv_image_width"), 1.0f / level_width);
311         glProgramUniform1f(motion_search_program, glGetUniformLocation(motion_search_program, "inv_image_height"), 1.0f / level_height);
312
313         GLuint flow_fbo;  // TODO: cleanup
314         glCreateFramebuffers(1, &flow_fbo);
315         glNamedFramebufferTexture(flow_fbo, GL_COLOR_ATTACHMENT0, flow_out_tex, 0);
316
317         glViewport(0, 0, width_patches, height_patches);
318         glBindFramebuffer(GL_FRAMEBUFFER, flow_fbo);
319         glBindVertexArray(motion_search_vao);
320         glUseProgram(motion_search_program);
321         glDrawArrays(GL_TRIANGLE_STRIP, 0, 4);
322 }
323
324 // Do “densification”, ie., upsampling of the flow patches to the flow field
325 // (the same size as the image at this level). We draw one quad per patch
326 // over its entire covered area (using instancing in the vertex shader),
327 // and then weight the contributions in the pixel shader by post-warp difference.
328 // This is equation (3) in the paper.
329 //
330 // We accumulate the flow vectors in the R/G channels (for u/v) and the total
331 // weight in the B channel. Dividing R and G by B gives the normalized values.
332 class Densify {
333 public:
334         Densify();
335         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);
336
337 private:
338         GLuint densify_vs_obj;
339         GLuint densify_fs_obj;
340         GLuint densify_program;
341         GLuint densify_vao;
342 };
343
344 Densify::Densify()
345 {
346         densify_vs_obj = compile_shader(read_file("densify.vert"), GL_VERTEX_SHADER);
347         densify_fs_obj = compile_shader(read_file("densify.frag"), GL_FRAGMENT_SHADER);
348         densify_program = link_program(densify_vs_obj, densify_fs_obj);
349
350         // Set up the VAO containing all the required position/texcoord data.
351         glCreateVertexArrays(1, &densify_vao);
352         glBindVertexArray(densify_vao);
353         glBindBuffer(GL_ARRAY_BUFFER, vertex_vbo);
354
355         GLint position_attrib = glGetAttribLocation(densify_program, "position");
356         glEnableVertexArrayAttrib(densify_vao, position_attrib);
357         glVertexAttribPointer(position_attrib, 2, GL_FLOAT, GL_FALSE, 0, BUFFER_OFFSET(0));
358 }
359
360 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)
361 {
362         glUseProgram(densify_program);
363
364         bind_sampler(densify_program, "image0_tex", 0, tex0_view, nearest_sampler);
365         bind_sampler(densify_program, "image1_tex", 1, tex1_view, linear_sampler);
366         bind_sampler(densify_program, "flow_tex", 2, flow_tex, nearest_sampler);
367
368         glProgramUniform1i(densify_program, glGetUniformLocation(densify_program, "width_patches"), width_patches);
369         glProgramUniform2f(densify_program, glGetUniformLocation(densify_program, "patch_size"),
370                 float(patch_size_pixels) / level_width,
371                 float(patch_size_pixels) / level_height);
372
373         float patch_spacing_x = float(level_width - patch_size_pixels) / (width_patches - 1);
374         float patch_spacing_y = float(level_height - patch_size_pixels) / (height_patches - 1);
375         glProgramUniform2f(densify_program, glGetUniformLocation(densify_program, "patch_spacing"),
376                 patch_spacing_x / level_width,
377                 patch_spacing_y / level_height);
378
379         GLuint dense_flow_fbo;  // TODO: cleanup
380         glCreateFramebuffers(1, &dense_flow_fbo);
381         glNamedFramebufferTexture(dense_flow_fbo, GL_COLOR_ATTACHMENT0, dense_flow_tex, 0);
382
383         glViewport(0, 0, level_width, level_height);
384         glEnable(GL_BLEND);
385         glBlendFunc(GL_ONE, GL_ONE);
386         glBindVertexArray(densify_vao);
387         glBindFramebuffer(GL_FRAMEBUFFER, dense_flow_fbo);
388         glDrawArraysInstanced(GL_TRIANGLE_STRIP, 0, 4, width_patches * height_patches);
389 }
390
391 int main(void)
392 {
393         if (SDL_Init(SDL_INIT_EVERYTHING) == -1) {
394                 fprintf(stderr, "SDL_Init failed: %s\n", SDL_GetError());
395                 exit(1);
396         }
397         SDL_GL_SetAttribute(SDL_GL_ALPHA_SIZE, 8);
398         SDL_GL_SetAttribute(SDL_GL_DEPTH_SIZE, 0);
399         SDL_GL_SetAttribute(SDL_GL_STENCIL_SIZE, 0);
400         SDL_GL_SetAttribute(SDL_GL_DOUBLEBUFFER, 1);
401
402         SDL_GL_SetAttribute(SDL_GL_CONTEXT_PROFILE_MASK, SDL_GL_CONTEXT_PROFILE_CORE);
403         SDL_GL_SetAttribute(SDL_GL_CONTEXT_MAJOR_VERSION, 4);
404         SDL_GL_SetAttribute(SDL_GL_CONTEXT_MINOR_VERSION, 5);
405         // SDL_GL_SetAttribute(SDL_GL_CONTEXT_FLAGS, SDL_GL_CONTEXT_DEBUG_FLAG);
406         SDL_Window *window = SDL_CreateWindow("OpenGL window",
407                         SDL_WINDOWPOS_UNDEFINED,
408                         SDL_WINDOWPOS_UNDEFINED,
409                         64, 64,
410                         SDL_WINDOW_OPENGL);
411         SDL_GLContext context = SDL_GL_CreateContext(window);
412         assert(context != nullptr);
413
414         // Load pictures.
415         GLuint tex0 = load_texture("test1499.pgm", WIDTH, HEIGHT);
416         GLuint tex1 = load_texture("test1500.pgm", WIDTH, HEIGHT);
417
418         // Make some samplers.
419         glCreateSamplers(1, &nearest_sampler);
420         glSamplerParameteri(nearest_sampler, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
421         glSamplerParameteri(nearest_sampler, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
422         glSamplerParameteri(nearest_sampler, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
423         glSamplerParameteri(nearest_sampler, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
424
425         glCreateSamplers(1, &linear_sampler);
426         glSamplerParameteri(linear_sampler, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
427         glSamplerParameteri(linear_sampler, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
428         glSamplerParameteri(linear_sampler, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
429         glSamplerParameteri(linear_sampler, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
430
431         glCreateSamplers(1, &mipmap_sampler);
432         glSamplerParameteri(mipmap_sampler, GL_TEXTURE_MIN_FILTER, GL_LINEAR_MIPMAP_NEAREST);
433         glSamplerParameteri(mipmap_sampler, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
434         glSamplerParameteri(mipmap_sampler, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
435         glSamplerParameteri(mipmap_sampler, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
436
437         float vertices[] = {
438                 0.0f, 1.0f,
439                 0.0f, 0.0f,
440                 1.0f, 1.0f,
441                 1.0f, 0.0f,
442         };
443         glCreateBuffers(1, &vertex_vbo);
444         glNamedBufferData(vertex_vbo, sizeof(vertices), vertices, GL_STATIC_DRAW);
445         glBindBuffer(GL_ARRAY_BUFFER, vertex_vbo);
446
447         // Initial flow is zero, 1x1.
448         GLuint initial_flow_tex;
449         glCreateTextures(GL_TEXTURE_2D, 1, &initial_flow_tex);
450         glTextureStorage2D(initial_flow_tex, 1, GL_RGB32F, 1, 1);
451
452         GLuint prev_level_flow_tex = initial_flow_tex;
453
454         Sobel sobel;
455         MotionSearch motion_search;
456         Densify densify;
457
458         for (int level = coarsest_level; level >= int(finest_level); --level) {
459                 int level_width = WIDTH >> level;
460                 int level_height = HEIGHT >> level;
461                 float patch_spacing_pixels = patch_size_pixels * (1.0f - patch_overlap_ratio);
462                 int width_patches = 1 + lrintf((level_width - patch_size_pixels) / patch_spacing_pixels);
463                 int height_patches = 1 + lrintf((level_height - patch_size_pixels) / patch_spacing_pixels);
464
465                 // Make sure we always read from the correct level; the chosen
466                 // mipmapping could otherwise be rather unpredictable, especially
467                 // during motion search.
468                 // TODO: create these beforehand, and stop leaking them.
469                 GLuint tex0_view, tex1_view;
470                 glGenTextures(1, &tex0_view);
471                 glTextureView(tex0_view, GL_TEXTURE_2D, tex0, GL_R8, level, 1, 0, 1);
472                 glGenTextures(1, &tex1_view);
473                 glTextureView(tex1_view, GL_TEXTURE_2D, tex1, GL_R8, level, 1, 0, 1);
474
475                 // Create a new texture; we could be fancy and render use a multi-level
476                 // texture, but meh.
477                 GLuint grad0_tex;
478                 glCreateTextures(GL_TEXTURE_2D, 1, &grad0_tex);
479                 glTextureStorage2D(grad0_tex, 1, GL_RG16F, level_width, level_height);
480
481                 // Find the derivative.
482                 sobel.exec(tex0_view, grad0_tex, level_width, level_height);
483
484                 // Motion search to find the initial flow. We use the flow from the previous
485                 // level (sampled bilinearly; no fancy tricks) as a guide, then search from there.
486
487                 // Create an output flow texture.
488                 GLuint flow_out_tex;
489                 glCreateTextures(GL_TEXTURE_2D, 1, &flow_out_tex);
490                 glTextureStorage2D(flow_out_tex, 1, GL_RG16F, width_patches, height_patches);
491
492                 // And draw.
493                 motion_search.exec(tex0_view, tex1_view, grad0_tex, prev_level_flow_tex, flow_out_tex, level_width, level_height, width_patches, height_patches);
494
495                 // Densification.
496
497                 // Set up an output texture.
498                 GLuint dense_flow_tex;
499                 glCreateTextures(GL_TEXTURE_2D, 1, &dense_flow_tex);
500                 //glTextureStorage2D(dense_flow_tex, 1, GL_RGB16F, level_width, level_height);
501                 glTextureStorage2D(dense_flow_tex, 1, GL_RGBA32F, level_width, level_height);
502
503                 // And draw.
504                 densify.exec(tex0_view, tex1_view, flow_out_tex, dense_flow_tex, level_width, level_height, width_patches, height_patches);
505
506                 // TODO: Variational refinement.
507
508                 prev_level_flow_tex = dense_flow_tex;
509         }
510
511         int level_width = WIDTH >> finest_level;
512         int level_height = HEIGHT >> finest_level;
513         unique_ptr<float[]> dense_flow(new float[level_width * level_height * 3]);
514         glGetTextureImage(prev_level_flow_tex, 0, GL_RGB, GL_FLOAT, level_width * level_height * 3 * sizeof(float), dense_flow.get());
515
516         FILE *fp = fopen("flow.ppm", "wb");
517         fprintf(fp, "P6\n%d %d\n255\n", level_width, level_height);
518         for (unsigned y = 0; y < unsigned(level_height); ++y) {
519                 int yy = level_height - y - 1;
520                 for (unsigned x = 0; x < unsigned(level_width); ++x) {
521                         float du = dense_flow[(yy * level_width + x) * 3 + 0];
522                         float dv = dense_flow[(yy * level_width + x) * 3 + 1];
523                         float w = dense_flow[(yy * level_width + x) * 3 + 2];
524
525                         du /= w;
526                         dv /= w;
527
528                         float angle = atan2(dv * level_width, du * level_height);
529                         float magnitude = min(hypot(du * level_width, dv * level_height) / 20.0, 1.0);
530                         
531                         // HSV to RGB (from Wikipedia). Saturation is 1.
532                         float c = magnitude;
533                         float h = (angle + M_PI) * 6.0 / (2.0 * M_PI);
534                         float X = c * (1.0 - fabs(fmod(h, 2.0) - 1.0));
535                         float r = 0.0f, g = 0.0f, b = 0.0f;
536                         if (h < 1.0f) {
537                                 r = c; g = X;
538                         } else if (h < 2.0f) {
539                                 r = X; g = c;
540                         } else if (h < 3.0f) {
541                                 g = c; b = X;
542                         } else if (h < 4.0f) {
543                                 g = X; b = c;
544                         } else if (h < 5.0f) {
545                                 r = X; b = c;
546                         } else if (h < 6.0f) {
547                                 r = c; b = X;
548                         } else {
549                                 // h is NaN, so black is fine.
550                         }
551                         float m = magnitude - c;
552                         r += m; g += m; b += m;
553                         r = max(min(r, 1.0f), 0.0f);
554                         g = max(min(g, 1.0f), 0.0f);
555                         b = max(min(b, 1.0f), 0.0f);
556                         putc(lrintf(r * 255.0f), fp);
557                         putc(lrintf(g * 255.0f), fp);
558                         putc(lrintf(b * 255.0f), fp);
559                 }
560         }
561         fclose(fp);
562
563         fprintf(stderr, "err = %d\n", glGetError());
564 }