Fix some errors in the smoothness term.

diff --git a/equations.frag b/equations.frag
index 6171a3729d6f862490e288e5125152227d6756b1..b8dbdde4216bb58e5381c7f0ed988b6af530e95b 100644 (file)
--- a/equations.frag
+++ b/equations.frag
@@ -4,17 +4,13 @@ in vec2 tc;
 out uvec4 equation;
 
 uniform sampler2D I_x_y_tex, I_t_tex;
-uniform sampler2D diff_flow_tex, flow_tex;
+uniform sampler2D diff_flow_tex, base_flow_tex;
 uniform sampler2D beta_0_tex;
 uniform sampler2D smoothness_x_tex, smoothness_y_tex;
 
 // TODO: Consider a specialized version for the case where we know that du = dv = 0,
 // since we run so few iterations.
 
-// This must be a macro, since the offset needs to be a constant expression.
-#define get_flow(x_offs, y_offs) \
-       (textureOffset(flow_tex, tc, ivec2((x_offs), (y_offs))).xy + \
-       textureOffset(diff_flow_tex, tc, ivec2((x_offs), (y_offs))).xy)
 
 void main()
 {
@@ -93,23 +89,18 @@ void main()
        float smooth_r = texture(smoothness_x_tex, tc).x;
        float smooth_d = textureOffset(smoothness_y_tex, tc, ivec2( 0, -1)).x;
        float smooth_u = texture(smoothness_y_tex, tc).x;
-       A11 -= smooth_l + smooth_r + smooth_d + smooth_u;
-       A22 -= smooth_l + smooth_r + smooth_d + smooth_u;
+       A11 += smooth_l + smooth_r + smooth_d + smooth_u;
+       A22 += smooth_l + smooth_r + smooth_d + smooth_u;
 
-       // Laplacian of (u0 + du, v0 + dv), sans the central term.
+       // Laplacian of (u0, v0).
        vec2 laplacian =
-               smooth_l * get_flow(-1, 0) +
-               smooth_r * get_flow(1, 0) +
-               smooth_d * get_flow(0, -1) +
-               smooth_u * get_flow(0, 1);
-       b1 -= laplacian.x;
-       b2 -= laplacian.y;
-
-       // The central term of the Laplacian, for (u0, v0) only.
-       // (The central term for (du, dv) is what we are solving for.)
-       vec2 central = (smooth_l + smooth_r + smooth_d + smooth_u) * texture(flow_tex, tc).xy;
-       b1 += central.x;
-       b2 += central.y;
+               smooth_l * textureOffset(base_flow_tex, tc, ivec2(-1,  0)).xy +
+               smooth_r * textureOffset(base_flow_tex, tc, ivec2( 1,  0)).xy +
+               smooth_d * textureOffset(base_flow_tex, tc, ivec2( 0, -1)).xy +
+               smooth_u * textureOffset(base_flow_tex, tc, ivec2( 0,  1)).xy -
+               (smooth_l + smooth_r + smooth_d + smooth_u) * texture(base_flow_tex, tc).xy;
+       b1 += laplacian.x;
+       b2 += laplacian.y;
 
        // Encode the equation down into four uint32s.
        equation.x = floatBitsToUint(1.0 / A11);

diff --git a/flow.cpp b/flow.cpp
index 9264123c24f3f0c70a7795f165a019491c8ad90b..15161988eb77ddab0c9eaf070385763f40038584 100644 (file)
--- a/flow.cpp
+++ b/flow.cpp
@@ -628,7 +628,7 @@ private:
        GLuint equations_vao;
 
        GLuint uniform_I_x_y_tex, uniform_I_t_tex;
-       GLuint uniform_diff_flow_tex, uniform_flow_tex;
+       GLuint uniform_diff_flow_tex, uniform_base_flow_tex;
        GLuint uniform_beta_0_tex;
        GLuint uniform_smoothness_x_tex, uniform_smoothness_y_tex;
 };
@@ -651,20 +651,20 @@ SetupEquations::SetupEquations()
        uniform_I_x_y_tex = glGetUniformLocation(equations_program, "I_x_y_tex");
        uniform_I_t_tex = glGetUniformLocation(equations_program, "I_t_tex");
        uniform_diff_flow_tex = glGetUniformLocation(equations_program, "diff_flow_tex");
-       uniform_flow_tex = glGetUniformLocation(equations_program, "flow_tex");
+       uniform_base_flow_tex = glGetUniformLocation(equations_program, "base_flow_tex");
        uniform_beta_0_tex = glGetUniformLocation(equations_program, "beta_0_tex");
        uniform_smoothness_x_tex = glGetUniformLocation(equations_program, "smoothness_x_tex");
        uniform_smoothness_y_tex = glGetUniformLocation(equations_program, "smoothness_y_tex");
 }
 
-void SetupEquations::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)
+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)
 {
        glUseProgram(equations_program);
 
        bind_sampler(equations_program, uniform_I_x_y_tex, 0, I_x_y_tex, nearest_sampler);
        bind_sampler(equations_program, uniform_I_t_tex, 1, I_t_tex, nearest_sampler);
        bind_sampler(equations_program, uniform_diff_flow_tex, 2, diff_flow_tex, nearest_sampler);
-       bind_sampler(equations_program, uniform_flow_tex, 3, flow_tex, nearest_sampler);
+       bind_sampler(equations_program, uniform_base_flow_tex, 3, base_flow_tex, nearest_sampler);
        bind_sampler(equations_program, uniform_beta_0_tex, 4, beta_0_tex, nearest_sampler);
        bind_sampler(equations_program, uniform_smoothness_x_tex, 5, smoothness_x_tex, smoothness_sampler);
        bind_sampler(equations_program, uniform_smoothness_y_tex, 6, smoothness_y_tex, smoothness_sampler);

diff --git a/smoothness.frag b/smoothness.frag
index 4dc06f3699dc526d804e82c879c4a4d2809243e1..583231ed90b233b69489facaa4d4a93641c8cce4 100644 (file)
--- a/smoothness.frag
+++ b/smoothness.frag
@@ -13,7 +13,7 @@ uniform sampler2D flow_tex, diff_flow_tex;
 
 float diffusivity(float u_x, float u_y, float v_x, float v_y)
 {
-       return -inversesqrt(u_x * u_x + u_y * u_y + v_x * v_x + v_y * v_y + eps_sq);
+       return inversesqrt(u_x * u_x + u_y * u_y + v_x * v_x + v_y * v_y + eps_sq);
 }
 
 void main()

diff --git a/sor.frag b/sor.frag
index f6ac2ee5adaf0780c67eb23783de58b45e17b222..dc4a5b8af4fdb19910819ffb39a96ee06246163f 100644 (file)
--- a/sor.frag
+++ b/sor.frag
@@ -23,10 +23,10 @@ void main()
        float smooth_r = texture(smoothness_x_tex, tc).x;
        float smooth_d = textureOffset(smoothness_y_tex, tc, ivec2( 0, -1)).x;
        float smooth_u = texture(smoothness_y_tex, tc).x;
-       b -= smooth_l * textureOffset(diff_flow_tex, tc, ivec2(-1,  0)).xy;
-       b -= smooth_r * textureOffset(diff_flow_tex, tc, ivec2( 1,  0)).xy;
-       b -= smooth_d * textureOffset(diff_flow_tex, tc, ivec2( 0, -1)).xy;
-       b -= smooth_u * textureOffset(diff_flow_tex, tc, ivec2( 0,  1)).xy;
+       b += smooth_l * textureOffset(diff_flow_tex, tc, ivec2(-1,  0)).xy;
+       b += smooth_r * textureOffset(diff_flow_tex, tc, ivec2( 1,  0)).xy;
+       b += smooth_d * textureOffset(diff_flow_tex, tc, ivec2( 0, -1)).xy;
+       b += smooth_u * textureOffset(diff_flow_tex, tc, ivec2( 0,  1)).xy;
 
        const float omega = 1.6;
        diff_flow = texture(diff_flow_tex, tc).xy;

diff --git a/variational_refinement.txt b/variational_refinement.txt
index 2f8a341f2b0280797adbb6195d2cd21722abf6f7..b94f7981636e59c3600a7b643d80cad3bff03dd3 100644 (file)
--- a/variational_refinement.txt
+++ b/variational_refinement.txt
@@ -335,15 +335,15 @@ operator; this factor has to be discretized together with u before we can treat
 it as a constant. So instead, we'll define it as a function (called the
 _diffusivity_ at the given point):
 
-  g(x, y) = -1 / sqrt(u_x² + u_y² + v_x² + v_y² + ε²) = 0
+  g(x, y) = 1 / sqrt(u_x² + u_y² + v_x² + v_y² + ε²) = 0
 
 which gives us
 
-  ∂/∂x ( g(x, y) u_x ) + ∂/∂y ( g(x, y) u_y ) = 0
+  - ∂/∂x ( g(x, y) u_x ) - ∂/∂y ( g(x, y) u_y ) = 0
 
 We'll also have a similar equation for minimizing v, of course:
 
-  ∂/∂x ( g(x, y) v_x ) + ∂/∂y ( g(x, y) v_y ) = 0
+  - ∂/∂x ( g(x, y) v_x ) - ∂/∂y ( g(x, y) v_y ) = 0
 
 There's no normalization term β here, unlike the other terms; DeepFlow2
 adds one, but we're not including it here.
@@ -364,11 +364,11 @@ differently in the numerics anyway.)
 
 This gives us:
 
-  ∂/∂x ( g(x, y) (u0 + du)_x ) + ∂/∂y ( g(x, y) (u0 + du)_y ) = 0
+  - ∂/∂x ( g(x, y) (u0 + du)_x ) - ∂/∂y ( g(x, y) (u0 + du)_y ) = 0
 
 or
 
-  ∂/∂x ( g(x, y) du_x ) + ∂/∂y ( g(x, y) du_y ) = - ∂/∂x ( g(x, y) u0_x ) - ∂/∂y ( g(x, y) u0_y )
+  - ∂/∂x ( g(x, y) du_x ) - ∂/∂y ( g(x, y) du_y ) = ∂/∂x ( g(x, y) u0_x ) + ∂/∂y ( g(x, y) u0_y )
 
 where the right-hand side is effectively a constant for these purposes
 (although it still needs to be calculated anew for each iteration,
@@ -399,21 +399,21 @@ du' and dv' will be zero.)
 
 Now back to our equations, as we look at practical implementation:
 
-  ∂/∂x ( g(x, y) du_x ) + ∂/∂y ( g(x, y) du_y ) = - ∂/∂x ( g(x, y) u0_x ) - ∂/∂y ( g(x, y) u0_y )
-  ∂/∂x ( g(x, y) dv_x ) + ∂/∂y ( g(x, y) dv_y ) = - ∂/∂x ( g(x, y) v0_x ) - ∂/∂y ( g(x, y) v0_y )
+  - ∂/∂x ( g(x, y) du_x ) - ∂/∂y ( g(x, y) du_y ) = ∂/∂x ( g(x, y) u0_x ) + ∂/∂y ( g(x, y) u0_y )
+  - ∂/∂x ( g(x, y) dv_x ) - ∂/∂y ( g(x, y) dv_y ) = ∂/∂x ( g(x, y) v0_x ) + ∂/∂y ( g(x, y) v0_y )
 
-We can discretize the left-hand and right-hand side identically, so let's
-look only at
+We can discretize the left-hand and right-hand side identically (they differ
+only in signs and in variable), so let's look only at
 
-  ∂/∂x ( g(x, y) du_x ) + ∂/∂y ( g(x, y) du_y )
+  - ∂/∂x ( g(x, y) du_x ) - ∂/∂y ( g(x, y) du_y )
 
 [Brox05] equation (2.14) (which refers to a 1998 book, although I couldn't
 immediately find the equation in question in that book) discretizes this as
 
-    1/2 (g(x+1, y) + g(x, y)) (du(x+1, y) - du(x, y))
-  - 1/2 (g(x-1, y) + g(x, y)) (du(x, y) - du(x-1, y))
-  + 1/2 (g(x, y+1) + g(x, y)) (du(x, y+1) - du(x, y))
-  - 1/2 (g(x, y-1) + g(x, y)) (du(x, y) - du(x, y-1))
+  - 1/2 (g(x+1, y) + g(x, y)) (du(x+1, y) - du(x, y))
+  + 1/2 (g(x-1, y) + g(x, y)) (du(x, y) - du(x-1, y))
+  - 1/2 (g(x, y+1) + g(x, y)) (du(x, y+1) - du(x, y))
+  + 1/2 (g(x, y-1) + g(x, y)) (du(x, y) - du(x, y-1))
 
 It also mentions that it would be better to sample g at the half-way points,
 e.g. g(x+0.5, y), but that begs the question exactly how we'd do that, and
@@ -422,7 +422,7 @@ DeepFlow doesn't seem to care, so we stick with their version.
 Now we can finally let g use the values of the flow (note that this is the
 actual flow u and v, not du and dv!) from the previous iteration, as before:
 
-  g(x, y) = -1 / sqrt(u'_x² + u'_y² + v'_x² + v'_y² + ε²) = 0
+  g(x, y) = 1 / sqrt(u'_x² + u'_y² + v'_x² + v'_y² + ε²)
 
 The single derivatives in g(x) are approximated by standard central differences
 (see https://en.wikipedia.org/wiki/Finite_difference_coefficient), e.g.
@@ -444,18 +444,18 @@ Let's now define a smoothness constant between the neighbors (x,y) and (x1,y1):
 Collecting all the du(x, y) terms of the discretized equation above,
 ignoring the right-hand side, which is just a constant for us anyway:
 
-    s(x+1, y) (du(x+1, y) - du(x, y))
-  - s(x-1, y) (du(x, y) - du(x-1, y))
-  + s(x, y+1) (du(x, y+1) - du(x, y))
-  - s(x, y-1) (du(x, y) - du(x, y-1)) = -C
+  - s(x+1, y) (du(x+1, y) - du(x, y))
+  + s(x-1, y) (du(x, y) - du(x-1, y))
+  - s(x, y+1) (du(x, y+1) - du(x, y))
+  + s(x, y-1) (du(x, y) - du(x, y-1)) = C
 
-    s(x+1, y) du(x+1, y) - s(x+1, y) du(x, y)
-  - s(x-1, y) du(x, y) + s(x-1, y) du(x-1, y)
-  + s(x, y+1) du(x, y+1) - s(x, y+1) du(x, y)
-  - s(x, y-1) du(x, y) + s(x, y-1) du(x, y-1) = -C
+  - s(x+1, y) du(x+1, y) + s(x+1, y) du(x, y)
+  + s(x-1, y) du(x, y) - s(x-1, y) du(x-1, y)
+  - s(x, y+1) du(x, y+1) + s(x, y+1) du(x, y)
+  + s(x, y-1) du(x, y) - s(x, y-1) du(x, y-1) = C
 
-  - (s(x+1, y) + s(x-1, y) + s(x, y+1) + s(x, y-1)) du(x, y) =
-  - s(x+1, y) du(x+1, y) - s(x-1, y) du(x-1, y) - s(x, y+1) du(x, y+1) - s(x, y-1) du(x, y-1) - C
+  (s(x+1, y) + s(x-1, y) + s(x, y+1) + s(x, y-1)) du(x, y) =
+  s(x+1, y) du(x+1, y) + s(x-1, y) du(x-1, y) + s(x, y+1) du(x, y+1) + s(x, y-1) du(x, y-1) + C
 
 It is interesting to note that if s = 1 uniformly, which would be the case
 without our penalizer Ψ(a²), we would have the familiar discrete Laplacian,