From c214d6919f05020b0695bccb8c948e6229c51375 Mon Sep 17 00:00:00 2001 From: "Steinar H. Gunderson" Date: Sat, 21 Jul 2018 01:30:14 +0200 Subject: [PATCH 1/1] Fix some errors in the smoothness term. --- equations.frag | 31 +++++++++-------------- flow.cpp | 8 +++--- smoothness.frag | 2 +- sor.frag | 8 +++--- variational_refinement.txt | 50 +++++++++++++++++++------------------- 5 files changed, 45 insertions(+), 54 deletions(-) diff --git a/equations.frag b/equations.frag index 6171a37..b8dbdde 100644 --- 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 9264123..1516198 100644 --- 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 4dc06f3..583231e 100644 --- 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 f6ac2ee..dc4a5b8 100644 --- 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 2f8a341..b94f798 100644 --- 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, -- 2.39.2