Add an implementation of sharpening by FIR Wiener filters.
authorSteinar H. Gunderson <sgunderson@bigfoot.com>
Wed, 10 Oct 2012 18:02:27 +0000 (20:02 +0200)
committerSteinar H. Gunderson <sgunderson@bigfoot.com>
Wed, 10 Oct 2012 18:02:27 +0000 (20:02 +0200)
Makefile
deconvolution_sharpen_effect.cpp [new file with mode: 0644]
deconvolution_sharpen_effect.frag [new file with mode: 0644]
deconvolution_sharpen_effect.h [new file with mode: 0644]

index 028b855..c33d20a 100644 (file)
--- 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 (file)
index 0000000..aa77dd8
--- /dev/null
@@ -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 <math.h>
+#include <assert.h>
+#include <Eigen/Dense>
+#include <Eigen/Cholesky>
+
+#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<int>(ya_min, 0);
+                       ya_max = std::min<int>(ya_max, a.rows() - 1);
+                       xa_min = std::max<int>(xa_min, 0);
+                       xa_max = std::min<int>(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<int>(ya_min, 0);
+                       ya_max = std::min<int>(ya_max, a.rows() - 1);
+                       xa_min = std::max<int>(xa_min, 0);
+                       xa_max = std::min<int>(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<MatrixXf> 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 (file)
index 0000000..7b90ac7
--- /dev/null
@@ -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 (file)
index 0000000..bc0bbd4
--- /dev/null
@@ -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)