From: Steinar H. Gunderson Date: Wed, 10 Oct 2012 18:02:27 +0000 (+0200) Subject: Add an implementation of sharpening by FIR Wiener filters. X-Git-Tag: 1.0~327 X-Git-Url: https://git.sesse.net/?p=movit;a=commitdiff_plain;h=42c35394ef92bb5179fc4879cb55b866fd422d28 Add an implementation of sharpening by FIR Wiener filters. --- diff --git a/Makefile b/Makefile index 028b855..c33d20a 100644 --- a/Makefile +++ b/Makefile @@ -1,6 +1,6 @@ CC=gcc CXX=g++ -CXXFLAGS=-Wall -g +CXXFLAGS=-Wall -g $(shell pkg-config --cflags eigen3 ) LDFLAGS=-lSDL -lSDL_image -lGL -lrt # Core. @@ -24,6 +24,7 @@ OBJS += diffusion_effect.o OBJS += glow_effect.o OBJS += mix_effect.o OBJS += resize_effect.o +OBJS += deconvolution_sharpen_effect.o OBJS += sandbox_effect.o test: $(OBJS) diff --git a/deconvolution_sharpen_effect.cpp b/deconvolution_sharpen_effect.cpp new file mode 100644 index 0000000..aa77dd8 --- /dev/null +++ b/deconvolution_sharpen_effect.cpp @@ -0,0 +1,438 @@ +// NOTE: Throughout, we use the symbol ⊙ for convolution. +// Since all of our signals are symmetrical, discrete correlation and convolution +// is the same operation, and so we won't make a difference in notation. + + +#include +#include +#include +#include + +#include "deconvolution_sharpen_effect.h" +#include "util.h" +#include "opengl.h" + +using namespace Eigen; + +DeconvolutionSharpenEffect::DeconvolutionSharpenEffect() + : R(5), + circle_radius(2.0f), + gaussian_radius(0.0f), + correlation(0.95f), + noise(0.01f) +{ + register_int("matrix_size", &R); + register_float("circle_radius", &circle_radius); + register_float("gaussian_radius", &gaussian_radius); + register_float("correlation", &correlation); + register_float("noise", &noise); +} + +std::string DeconvolutionSharpenEffect::output_fragment_shader() +{ + char buf[256]; + sprintf(buf, "#define R %u\n", R); + return buf + read_file("deconvolution_sharpen_effect.frag"); +} + +namespace { + +// Integral of sqrt(r² - x²) dx over x=0..a. +float circle_integral(float a, float r) +{ + assert(a >= 0.0f); + if (a <= 0.0f) { + return 0.0f; + } + if (a >= r) { + return 0.25f * M_PI * r * r; + } + return 0.5f * (a * sqrt(r*r - a*a) + r*r * asin(a / r)); +} + +// Yields the impulse response of a circular blur with radius r. +// We basically look at each element as a square centered around (x,y), +// and figure out how much of its area is covered by the circle. +float circle_impulse_response(int x, int y, float r) +{ + if (r < 1e-3) { + // Degenerate case: radius = 0 yields the impulse response. + return (x == 0 && y == 0) ? 1.0f : 0.0f; + } + + // Find the extents of this cell. Due to symmetry, we can cheat a bit + // and pretend we're always in the upper-right quadrant, except when + // we're right at an axis crossing (x = 0 or y = 0), in which case we + // simply use the evenness of the function; shrink the cell, make + // the calculation, and down below we'll normalize by the cell's area. + float min_x, max_x, min_y, max_y; + if (x == 0) { + min_x = 0.0f; + max_x = 0.5f; + } else { + min_x = abs(x) - 0.5f; + max_x = abs(x) + 0.5f; + } + if (y == 0) { + min_y = 0.0f; + max_y = 0.5f; + } else { + min_y = abs(y) - 0.5f; + max_y = abs(y) + 0.5f; + } + assert(min_x >= 0.0f && max_x >= 0.0f); + assert(min_y >= 0.0f && max_y >= 0.0f); + + float cell_height = max_y - min_y; + float cell_width = max_x - min_x; + + if (min_x * min_x + min_y * min_y > r * r) { + // Lower-left corner is outside the circle, so the entire cell is. + return 0.0f; + } + if (max_x * max_x + max_y * max_y < r * r) { + // Upper-right corner is inside the circle, so the entire cell is. + return 1.0f; + } + + // OK, so now we know the cell is partially covered by the circle: + // + // \ . + // ------------- + // |####|#\ | + // |####|##| | + // ------------- + // A ###| + // ###| + // + // The edge of the circle is defined by x² + y² = r², + // or x = sqrt(r² - y²) (since x is nonnegative). + // Find out where the curve crosses our given y values. + float mid_x1 = (max_y >= r) ? min_x : sqrt(r * r - max_y * max_y); + float mid_x2 = sqrt(r * r - min_y * min_y); + if (mid_x1 < min_x) { + mid_x1 = min_x; + } + if (mid_x2 > max_x) { + mid_x2 = max_x; + } + assert(mid_x1 >= min_x); + assert(mid_x2 >= mid_x1); + assert(max_x >= mid_x2); + + // The area marked A in the figure above. + float covered_area = cell_height * (mid_x1 - min_x); + + // The area marked B in the figure above. Note that the integral gives the entire + // shaded space down to zero, so we need to subtract the rectangle that does not + // belong to our cell. + covered_area += circle_integral(mid_x2, r) - circle_integral(mid_x1, r); + covered_area -= min_y * (mid_x2 - mid_x1); + + assert(covered_area <= cell_width * cell_height); + return covered_area / (cell_width * cell_height); +} + +// Compute a ⊙ b. Note that we compute the “full” convolution, +// ie., our matrix will be big enough to hold every nonzero element of the result. +MatrixXf convolve(const MatrixXf &a, const MatrixXf &b) +{ + MatrixXf result(a.rows() + b.rows() - 1, a.cols() + b.cols() - 1); + for (int yr = 0; yr < result.rows(); ++yr) { + for (int xr = 0; xr < result.cols(); ++xr) { + float sum = 0.0f; + + // Given that x_b = x_r - x_a, find the values of x_a where + // x_a is in [0, a_cols> and x_b is in [0, b_cols>. (y is similar.) + // + // The second demand gives: + // + // 0 <= x_r - x_a < b_cols + // 0 >= x_a - x_r > -b_cols + // x_r >= x_a > x_r - b_cols + int ya_min = yr - b.rows() + 1; + int ya_max = yr; + int xa_min = xr - b.rows() + 1; + int xa_max = xr; + + // Now fit to the first demand. + ya_min = std::max(ya_min, 0); + ya_max = std::min(ya_max, a.rows() - 1); + xa_min = std::max(xa_min, 0); + xa_max = std::min(xa_max, a.cols() - 1); + + assert(ya_max >= ya_min); + assert(xa_max >= xa_min); + + for (int ya = ya_min; ya <= ya_max; ++ya) { + for (int xa = xa_min; xa <= xa_max; ++xa) { + sum += a(ya, xa) * b(yr - ya, xr - xa); + } + } + + result(yr, xr) = sum; + } + } + return result; +} + +// Similar to convolve(), but instead of assuming every element outside +// of b is zero, we make no such assumption and instead return only the +// elements where we know the right answer. (This is the only difference +// between the two.) +// This is the same as conv2(a, b, 'valid') in Octave. +// +// a must be the larger matrix of the two. +MatrixXf central_convolve(const MatrixXf &a, const MatrixXf &b) +{ + assert(a.rows() >= b.rows()); + assert(a.cols() >= b.cols()); + MatrixXf result(a.rows() - b.rows() + 1, a.cols() - b.cols() + 1); + for (int yr = b.rows() - 1; yr < result.rows() + b.rows() - 1; ++yr) { + for (int xr = b.cols() - 1; xr < result.cols() + b.cols() - 1; ++xr) { + float sum = 0.0f; + + // Given that x_b = x_r - x_a, find the values of x_a where + // x_a is in [0, a_cols> and x_b is in [0, b_cols>. (y is similar.) + // + // The second demand gives: + // + // 0 <= x_r - x_a < b_cols + // 0 >= x_a - x_r > -b_cols + // x_r >= x_a > x_r - b_cols + int ya_min = yr - b.rows() + 1; + int ya_max = yr; + int xa_min = xr - b.rows() + 1; + int xa_max = xr; + + // Now fit to the first demand. + ya_min = std::max(ya_min, 0); + ya_max = std::min(ya_max, a.rows() - 1); + xa_min = std::max(xa_min, 0); + xa_max = std::min(xa_max, a.cols() - 1); + + assert(ya_max >= ya_min); + assert(xa_max >= xa_min); + + for (int ya = ya_min; ya <= ya_max; ++ya) { + for (int xa = xa_min; xa <= xa_max; ++xa) { + sum += a(ya, xa) * b(yr - ya, xr - xa); + } + } + + result(yr - b.rows() + 1, xr - b.cols() + 1) = sum; + } + } + return result; +} + +void print_matrix(const MatrixXf &m) +{ + for (int y = 0; y < m.rows(); ++y) { + for (int x = 0; x < m.cols(); ++x) { + printf("%7.4f ", m(x, y)); + } + printf("\n"); + } +} + +} // namespace + +void DeconvolutionSharpenEffect::set_gl_state(GLuint glsl_program_num, const std::string &prefix, unsigned *sampler_num) +{ + Effect::set_gl_state(glsl_program_num, prefix, sampler_num); + + assert(R >= 1); + assert(R <= 25); // Same limit as Refocus. + + printf("circular blur radius: %5.3f\n", circle_radius); + printf("gaussian blur radius: %5.3f\n", gaussian_radius); + printf("correlation: %5.3f\n", correlation); + printf("noise factor: %5.3f\n", noise); + printf("\n"); + + // Figure out the impulse response for the circular part of the blur. + MatrixXf circ_h(2 * R + 1, 2 * R + 1); + for (int y = -R; y <= R; ++y) { + for (int x = -R; x <= R; ++x) { + circ_h(y + R, x + R) = circle_impulse_response(x, y, circle_radius); + } + } + + // Same, for the Gaussian part of the blur. We make this a lot larger + // since we're going to convolve with it soon, and it has infinite support + // (see comments for central_convolve()). + MatrixXf gaussian_h(4 * R + 1, 4 * R + 1); + for (int y = -2 * R; y <= 2 * R; ++y) { + for (int x = -2 * R; x <= 2 * R; ++x) { + float val; + if (gaussian_radius < 1e-3) { + val = (x == 0 && y == 0) ? 1.0f : 0.0f; + } else { + float z = hypot(x, y) / gaussian_radius; + val = exp(-z * z); + } + gaussian_h(y + 2 * R, x + 2 * R) = val; + } + } + + // h, the (assumed) impulse response that we're trying to invert. + MatrixXf h = central_convolve(gaussian_h, circ_h); + assert(h.rows() == 2 * R + 1); + assert(h.cols() == 2 * R + 1); + + // Normalize the impulse response. + float sum = 0.0f; + for (int y = 0; y < 2 * R + 1; ++y) { + for (int x = 0; x < 2 * R + 1; ++x) { + sum += h(y, x); + } + } + for (int y = 0; y < 2 * R + 1; ++y) { + for (int x = 0; x < 2 * R + 1; ++x) { + h(y, x) /= sum; + } + } + + // r_uu, the (estimated/assumed) autocorrelation of the input signal (u). + // The signal is modelled a standard autoregressive process with the + // given correlation coefficient. + // + // We have to take a bit of care with the size of this matrix. + // The pow() function naturally has an infinite support (except for the + // degenerate case of correlation=0), but we have to chop it off + // somewhere. Since we convolve it with a 4*R+1 large matrix below, + // we need to make it twice as big as that, so that we have enough + // data to make r_vv valid. (central_convolve() effectively enforces + // that we get at least the right size.) + MatrixXf r_uu(8 * R + 1, 8 * R + 1); + for (int y = -4 * R; y <= 4 * R; ++y) { + for (int x = -4 * R; x <= 4 * R; ++x) { + r_uu(x + 4 * R, y + 4 * R) = pow(correlation, hypot(x, y)); + } + } + + // Estimate r_vv, the autocorrelation of the output signal v. + // Since we know that v = h ⊙ u and both are symmetrical, + // convolution and correlation are the same, and + // r_vv = v ⊙ v = (h ⊙ u) ⊙ (h ⊙ u) = (h ⊙ h) ⊙ r_uu. + MatrixXf r_vv = central_convolve(r_uu, convolve(h, h)); + assert(r_vv.rows() == 4 * R + 1); + assert(r_vv.cols() == 4 * R + 1); + + // Similarly, r_uv = u ⊙ v = u ⊙ (h ⊙ u) = h ⊙ r_uu. + //MatrixXf r_uv = central_convolve(r_uu, h).block(2 * R, 2 * R, 2 * R + 1, 2 * R + 1); + MatrixXf r_uu_center = r_uu.block(2 * R, 2 * R, 4 * R + 1, 4 * R + 1); + MatrixXf r_uv = central_convolve(r_uu_center, h); + assert(r_uv.rows() == 2 * R + 1); + assert(r_uv.cols() == 2 * R + 1); + + // Add the noise term (we assume the noise is uncorrelated, + // so it only affects the central element). + r_vv(2 * R, 2 * R) += noise; + + // Now solve the Wiener-Hopf equations to find the deconvolution kernel g. + // Most texts show this only for the simpler 1D case: + // + // [ r_vv(0) r_vv(1) r_vv(2) ... ] [ g(0) ] [ r_uv(0) ] + // [ r_vv(-1) r_vv(0) ... ] [ g(1) ] = [ r_uv(1) ] + // [ r_vv(-2) ... ] [ g(2) ] [ r_uv(2) ] + // [ ... ] [ g(3) ] [ r_uv(3) ] + // + // (Since r_vv is symmetrical, we can drop the minus signs.) + // + // Generally, row i of the matrix contains (dropping _vv for brevity): + // + // [ r(0-i) r(1-i) r(2-i) ... ] + // + // However, we have the 2D case. We flatten the vectors out to + // 1D quantities; this means we must think of the row number + // as a pair instead of as a scalar. Row (i,j) then contains: + // + // [ r(0-i,0-j) r(1-i,0-j) r(2-i,0-j) ... r(0-i,1-j) r_(1-i,1-j) r(2-i,1-j) ... ] + // + // g and r_uv are flattened in the same fashion. + // + // Note that even though this matrix is block Toeplitz, it is _not_ Toeplitz, + // and thus can not be inverted through the standard Levinson-Durbin method. + // There exists a block Levinson-Durbin method, which we may or may not + // want to use later. (Eigen's solvers are fast enough that for big matrices, + // the convolution operation and not the matrix solving is the bottleneck.) + // + // One thing we definitely want to use, though, is the symmetry properties. + // Since we know that g(i, j) = g(|i|, |j|), we can reduce the amount of + // unknowns to about 1/4th of the total size. The method is quite simple, + // as can be seen from the following toy equation system: + // + // A x0 + B x1 + C x2 = y0 + // D x0 + E x1 + F x2 = y1 + // G x0 + H x1 + I x2 = y2 + // + // If we now know that e.g. x0=x1 and y0=y1, we can rewrite this to + // + // (A+B+D+E) x0 + (C+F) x2 = 2 y0 + // (G+H) x0 + I x2 = y2 + // + // This both increases accuracy and provides us with a very nice speed + // boost. We could have gone even further and went for 8-way symmetry + // like the shader does, but this is good enough right now. + MatrixXf M(MatrixXf::Zero((R + 1) * (R + 1), (R + 1) * (R + 1))); + MatrixXf r_uv_flattened(MatrixXf::Zero((R + 1) * (R + 1), 1)); + for (int outer_i = 0; outer_i < 2 * R + 1; ++outer_i) { + int folded_outer_i = abs(outer_i - R); + for (int outer_j = 0; outer_j < 2 * R + 1; ++outer_j) { + int folded_outer_j = abs(outer_j - R); + int row = folded_outer_i * (R + 1) + folded_outer_j; + for (int inner_i = 0; inner_i < 2 * R + 1; ++inner_i) { + int folded_inner_i = abs(inner_i - R); + for (int inner_j = 0; inner_j < 2 * R + 1; ++inner_j) { + int folded_inner_j = abs(inner_j - R); + int col = folded_inner_i * (R + 1) + folded_inner_j; + M(row, col) += r_vv((inner_i - R) - (outer_i - R) + 2 * R, + (inner_j - R) - (outer_j - R) + 2 * R); + } + } + r_uv_flattened(row) += r_uv(outer_i, outer_j); + } + } + + LLT llt(M); + MatrixXf g_flattened = llt.solve(r_uv_flattened); + assert(g_flattened.rows() == (R + 1) * (R + 1)), + assert(g_flattened.cols() == 1); + + // Normalize and de-flatten the deconvolution matrix. + MatrixXf g(R + 1, R + 1); + sum = 0.0f; + for (int i = 0; i < g_flattened.rows(); ++i) { + int y = i / (R + 1); + int x = i % (R + 1); + if (y == 0 && x == 0) { + sum += g_flattened(i); + } else if (y == 0 || x == 0) { + sum += 2.0f * g_flattened(i); + } else { + sum += 4.0f * g_flattened(i); + } + } + for (int i = 0; i < g_flattened.rows(); ++i) { + int y = i / (R + 1); + int x = i % (R + 1); + g(y, x) = g_flattened(i) / sum; + } + + // Now encode it as uniforms, and pass it on to the shader. + // (Actually the shader only uses about half of the elements.) + float samples[4 * (R + 1) * (R + 1)]; + for (int y = 0; y <= R; ++y) { + for (int x = 0; x <= R; ++x) { + int i = y * (R + 1) + x; + samples[i * 4 + 0] = x / float(width); + samples[i * 4 + 1] = y / float(height); + samples[i * 4 + 2] = g(y, x); + samples[i * 4 + 3] = 0.0f; + } + } + + set_uniform_vec4_array(glsl_program_num, prefix, "samples", samples, R * R); +} diff --git a/deconvolution_sharpen_effect.frag b/deconvolution_sharpen_effect.frag new file mode 100644 index 0000000..7b90ac7 --- /dev/null +++ b/deconvolution_sharpen_effect.frag @@ -0,0 +1,66 @@ +uniform vec4 PREFIX(samples)[(R + 1) * (R + 1)]; + +vec4 FUNCNAME(vec2 tc) { + // The full matrix has five different symmetry cases, that look like this: + // + // D * * C * * D + // * D * C * D * + // * * D C D * * + // B B B A B B B + // * * D C D * * + // * D * C * D * + // D * * C * * D + // + // We only store the lower-right part of the matrix: + // + // A B B + // C D * + // C * D + + // Case A: Top-left sample has no symmetry. + vec4 sum = PREFIX(samples)[0].z * INPUT(tc); + + // Case B: Uppermost samples have left/right symmetry. + for (int x = 1; x <= R; ++x) { + vec4 sample = PREFIX(samples)[x]; + sum += sample.z * (INPUT(tc - sample.xy) + INPUT(tc + sample.xy)); + } + + // Case C: Leftmost samples have top/bottom symmetry. + for (int y = 1; y <= R; ++y) { + vec4 sample = PREFIX(samples)[y * (R + 1)]; + sum += sample.z * (INPUT(tc - sample.xy) + INPUT(tc + sample.xy)); + } + + // Case D: Diagonal samples have four-way symmetry. + for (int xy = 1; xy <= R; ++xy) { + vec4 sample = PREFIX(samples)[xy * (R + 1) + xy]; + + vec4 local_sum = INPUT(tc - sample.xy) + INPUT(tc + sample.xy); + sample.y = -sample.y; + local_sum += INPUT(tc - sample.xy) + INPUT(tc + sample.xy); + + sum += sample.z * local_sum; + } + + // Case *: All other samples have eight-way symmetry. + for (int y = 1; y <= R; ++y) { + for (int x = y + 1; x <= R; ++x) { + vec4 sample = PREFIX(samples)[y * (R + 1) + x]; + vec2 mirror_sample = vec2(sample.x, -sample.y); + + vec4 local_sum = INPUT(tc - sample.xy) + INPUT(tc + sample.xy); + local_sum += INPUT(tc - mirror_sample.xy) + INPUT(tc + mirror_sample.xy); + + sample.xy = sample.yx; + mirror_sample.xy = mirror_sample.yx; + + local_sum += INPUT(tc - sample.xy) + INPUT(tc + sample.xy); + local_sum += INPUT(tc - mirror_sample.xy) + INPUT(tc + mirror_sample.xy); + + sum += sample.z * local_sum; + } + } + + return sum; +} diff --git a/deconvolution_sharpen_effect.h b/deconvolution_sharpen_effect.h new file mode 100644 index 0000000..bc0bbd4 --- /dev/null +++ b/deconvolution_sharpen_effect.h @@ -0,0 +1,56 @@ +#ifndef _DECONVOLUTION_SHARPEN_EFFECT_H +#define _DECONVOLUTION_SHARPEN_EFFECT_H 1 + +// DeconvolutionSharpenEffect is an effect that sharpens by way of deconvolution +// (i.e., trying to reverse the blur kernel, as opposed to just boosting high +// frequencies), more specifically by FIR Wiener filters. It is the same +// algorithm as used by the (now largely abandoned) Refocus plug-in for GIMP, +// and I suspect the same as in Photoshop's “Smart Sharpen” filter. +// The implementation is, however, distinct from either. +// +// The effect gives generally better results than unsharp masking, but can be very +// GPU intensive, and requires a fair bit of tweaking to get good results without +// ringing and/or excessive noise. It should be mentioned that for the larger +// convolutions (e.g. R approaching 10), we should probably move to FFT-based +// convolution algorithms, especially as Mesa's shader compiler starts having +// problems compiling our shader. +// +// We follow the same book as Refocus was implemented from, namely +// +// Jain, Anil K.: “Fundamentals of Digital Image Processing”, Prentice Hall, 1988. + +#include "effect.h" + +class DeconvolutionSharpenEffect : public Effect { +public: + DeconvolutionSharpenEffect(); + virtual std::string effect_type_id() const { return "DeconvolutionSharpenEffect"; } + std::string output_fragment_shader(); + + virtual void inform_input_size(unsigned input_num, unsigned width, unsigned height) + { + this->width = width; + this->height = height; + } + + void set_gl_state(GLuint glsl_program_num, const std::string &prefix, unsigned *sampler_num); + +private: + // Input size. + unsigned width, height; + + // The maximum radius of the (de)convolution kernel. + // Note that since this extends both ways, and we also have a center element, + // the actual convolution matrix will be (2R + 1) x (2R + 1). + // + // Must match the definition in the shader, and as such, cannot be set once + // the chain has been finalized. + int R; + + // The parameters. Typical OK values are circle_radius = 2, gaussian_radius = 0 + // (ie., blur is assumed to be a 2px circle), correlation = 0.95, and noise = 0.01. + // Note that once the radius starts going too far past R, you will get nonsensical results. + float circle_radius, gaussian_radius, correlation, noise; +}; + +#endif // !defined(_DECONVOLUTION_SHARPEN_EFFECT_H)