Add the densification step, in theory now providing us with flow. But still too buggy...

diff --git a/densify.frag b/densify.frag
new file mode 100644 (file)
index 0000000..d16d630
--- /dev/null
+++ b/densify.frag
@@ -0,0 +1,21 @@
+#version 450 core
+
+in vec2 image_pos;
+flat in vec2 flow_du;
+out vec3 flow_contribution;
+
+uniform sampler2D image0_tex, image1_tex;
+
+void main()
+{
+       // Equation (3) from the paper. We're using additive blending, so the
+       // sum will happen automatically for us, and normalization happens on
+       // next read.
+       //
+       // Note that equation (2) says 1 for the minimum error, but the code says 2.0.
+       // And it says L2 norm, but really, the code does absolute value even for
+       // L2 error norm (it uses a square root formula for L1 norm).
+       float diff = texture(image0_tex, image_pos).x - texture(image1_tex, image_pos + flow_du).x;
+       float weight = 1.0 / max(abs(diff), 2.0 / 255.0);
+       flow_contribution = vec3(flow_du.x * weight, flow_du.y * weight, weight);
+}

diff --git a/densify.vert b/densify.vert
new file mode 100644 (file)
index 0000000..6a02041
--- /dev/null
+++ b/densify.vert
@@ -0,0 +1,32 @@
+#version 450 core
+
+in vec2 position;
+out vec2 image_pos;
+flat out vec2 flow_du;
+
+uniform int width_patches;
+uniform vec2 patch_size;  // In 0..1 coordinates.
+uniform vec2 patch_spacing;  // In 0..1 coordinates.
+uniform sampler2D flow_tex;
+
+void main()
+{
+       int patch_x = gl_InstanceID % width_patches;
+       int patch_y = gl_InstanceID / width_patches;
+
+       // TODO: Lock the bottom left of the patch to an integer number of pixels?
+
+       image_pos = patch_spacing * ivec2(patch_x, patch_y) + patch_size * position;
+
+       // Find the flow value for this patch, and send it on to the fragment shader.
+       vec2 patch_center = patch_spacing * ivec2(patch_x, patch_y) + patch_size * vec2(0.5, 0.5);
+       flow_du = texture(flow_tex, patch_center).xy;
+
+       // The result of glOrtho(0.0, 1.0, 0.0, 1.0, 0.0, 1.0) is:
+       //
+       //   2.000  0.000  0.000 -1.000
+       //   0.000  2.000  0.000 -1.000
+       //   0.000  0.000 -2.000 -1.000
+       //   0.000  0.000  0.000  1.000
+       gl_Position = vec4(2.0 * image_pos.x - 1.0, 2.0 * image_pos.y - 1.0, -1.0, 1.0);
+}

diff --git a/flow.cpp b/flow.cpp
index f4fbb74d67903c5ee1d422977b93636ae3026a79..7eec660b30c26e1b3b7783cce7b92dd528c77c3a 100644 (file)
--- a/flow.cpp
+++ b/flow.cpp
@@ -235,6 +235,10 @@ int main(void)
        GLuint sobel_fs_obj = compile_shader(read_file("sobel.frag"), GL_FRAGMENT_SHADER);
        GLuint sobel_program = link_program(sobel_vs_obj, sobel_fs_obj);
 
+       GLuint densify_vs_obj = compile_shader(read_file("densify.vert"), GL_VERTEX_SHADER);
+       GLuint densify_fs_obj = compile_shader(read_file("densify.frag"), GL_FRAGMENT_SHADER);
+       GLuint densify_program = link_program(densify_vs_obj, densify_fs_obj);
+
        // Make some samplers.
        GLuint nearest_sampler;
        glCreateSamplers(1, &nearest_sampler);
@@ -329,6 +333,7 @@ int main(void)
        glViewport(0, 0, level_width, level_height);
        glBindFramebuffer(GL_FRAMEBUFFER, grad0_fbo);
        glUseProgram(sobel_program);
+       glDisable(GL_BLEND);
        glDrawArrays(GL_TRIANGLE_STRIP, 0, 4);
 
        // Motion search to find the initial flow.
@@ -384,5 +389,111 @@ int main(void)
        glUseProgram(motion_search_program);
        glDrawArrays(GL_TRIANGLE_STRIP, 0, 4);
 
+       // Do “densification”, ie., upsampling of the flow patches to the flow field
+       // (the same size as the image at this level). We draw one quad per patch
+       // over its entire covered area (using instancing in the vertex shader),
+       // and then weight the contributions in the pixel shader by post-warp difference.
+       // This is equation (3) in the paper.
+       //
+       // We accumulate the flow vectors in the R/G channels (for u/v) and the total
+       // weight in the B channel. Dividing R and G by B gives the normalized values.
+
+       // Set up an output texture.
+       GLuint dense_flow_tex;
+       glCreateTextures(GL_TEXTURE_2D, 1, &dense_flow_tex);
+       //glTextureStorage2D(dense_flow_tex, 1, GL_RGB16F, level_width, level_height);
+       glTextureStorage2D(dense_flow_tex, 1, GL_RGBA32F, level_width, level_height);
+
+       GLuint dense_flow_fbo;
+       glCreateFramebuffers(1, &dense_flow_fbo);
+       glNamedFramebufferTexture(dense_flow_fbo, GL_COLOR_ATTACHMENT0, dense_flow_tex, 0);
+
+       glUseProgram(densify_program);
+
+       bind_sampler(densify_program, "image0_tex", 0, tex0_view, nearest_sampler);
+       bind_sampler(densify_program, "image1_tex", 1, tex1_view, linear_sampler);
+       bind_sampler(densify_program, "flow_tex", 2, flow_out_tex, nearest_sampler);
+
+       glProgramUniform1i(densify_program, glGetUniformLocation(densify_program, "width_patches"), width_patches);
+       glProgramUniform2f(densify_program, glGetUniformLocation(densify_program, "patch_size"),
+               float(patch_size_pixels) / level_width,
+               float(patch_size_pixels) / level_height);
+
+       float patch_spacing_x = float(level_width - patch_size_pixels) / (width_patches - 1);
+       float patch_spacing_y = float(level_height - patch_size_pixels) / (height_patches - 1);
+       glProgramUniform2f(densify_program, glGetUniformLocation(densify_program, "patch_spacing"),
+               patch_spacing_x / level_width,
+               patch_spacing_y / level_height);
+
+       // Set up the VAO containing all the required position/texcoord data.
+       GLuint densify_vao;
+        glCreateVertexArrays(1, &densify_vao);
+        glBindVertexArray(densify_vao);
+       glBindBuffer(GL_ARRAY_BUFFER, vertex_vbo);
+
+       position_attrib = glGetAttribLocation(densify_program, "position");
+       glEnableVertexArrayAttrib(densify_vao, position_attrib);
+       glVertexAttribPointer(position_attrib, 2, GL_FLOAT, GL_FALSE, 0, BUFFER_OFFSET(0));
+
+       glBindBuffer(GL_ARRAY_BUFFER, 0);
+
+       // And draw.
+       glViewport(0, 0, level_width, level_height);
+       glEnable(GL_BLEND);
+       glBlendFunc(GL_ONE, GL_ONE);
+       glBindFramebuffer(GL_FRAMEBUFFER, dense_flow_fbo);
+       glDrawArraysInstanced(GL_TRIANGLE_STRIP, 0, 4, width_patches * height_patches);
+
+       // TODO: Variational refinement.
+
+       unique_ptr<float[]> dense_flow(new float[level_width * level_height * 3]);
+       glGetTextureImage(dense_flow_tex, 0, GL_RGB, GL_FLOAT, level_width * level_height * 3 * sizeof(float), dense_flow.get());
+
+       FILE *fp = fopen("flow.ppm", "wb");
+       fprintf(fp, "P6\n%d %d\n255\n", level_width, level_height);
+       for (unsigned y = 0; y < level_height; ++y) {
+               for (unsigned x = 0; x < level_width; ++x) {
+                       float du = dense_flow[(y * level_width + x) * 3 + 0];
+                       float dv = dense_flow[(y * level_width + x) * 3 + 1];
+                       float w = dense_flow[(y * level_width + x) * 3 + 2];
+
+                       du /= w;
+                       dv /= w;
+
+                       float angle = atan2(dv, du);
+                       float magnitude = min(hypot(du * level_width, dv * level_height) / 20.0, 1.0);
+                       
+                       // HSV to RGB (from Wikipedia). Saturation is 1.
+                       float c = magnitude;
+                       float h = angle * 6.0 / (2.0 * M_PI);
+                       float X = c * (1.0 - (fmod(h, 2.0) - 1.0));
+                       float r = 0.0f, g = 0.0f, b = 0.0f;
+                       if (h < 1.0f) {
+                               r = c; g = X;
+                       } else if (h < 2.0f) {
+                               r = X; g = c;
+                       } else if (h < 3.0f) {
+                               g = c; b = X;
+                       } else if (h < 4.0f) {
+                               g = X; b = c;
+                       } else if (h < 5.0f) {
+                               r = X; b = c;
+                       } else if (h < 6.0f) {
+                               r = c; b = X;
+                       } else {
+                               // h is NaN, so black is fine.
+                       }
+                       float m = magnitude - c;
+                       r += m, g += m, b += m;
+                       r = max(min(r, 1.0f), 0.0f);
+                       g = max(min(g, 1.0f), 0.0f);
+                       b = max(min(b, 1.0f), 0.0f);
+                       putc(lrintf(r * 255.0f), fp);
+                       putc(lrintf(g * 255.0f), fp);
+                       putc(lrintf(b * 255.0f), fp);
+               }
+       }
+       fclose(fp);
+
        fprintf(stderr, "err = %d\n", glGetError());
 }

diff --git a/motion_search.frag b/motion_search.frag
index 616e2080dc1084342dd1723c283ba9deab755fc9..98ebff81b814357cb471e864bbd4c70baa0a5991 100644 (file)
--- a/motion_search.frag
+++ b/motion_search.frag
@@ -107,7 +107,10 @@ void main()
                u += H_inv * du * vec2(inv_image_width, inv_image_height);
        }
 
-       // TODO: reject if moving too far
+       // Reject if we moved too far.
+       if (length((u - initial_u) * vec2(image_width, image_height)) > patch_size) {
+               u = initial_u;
+       }
 
        out_flow = u;
 }