1 // NOTE: Throughout, we use the symbol ⊙ for convolution.
2 // Since all of our signals are symmetrical, discrete correlation and convolution
3 // is the same operation, and so we won't make a difference in notation.
6 #include <Eigen/Cholesky>
15 #include "deconvolution_sharpen_effect.h"
16 #include "effect_util.h"
19 using namespace Eigen;
24 DeconvolutionSharpenEffect::DeconvolutionSharpenEffect()
27 gaussian_radius(0.0f),
31 last_circle_radius(-1.0f),
32 last_gaussian_radius(-1.0f),
33 last_correlation(-1.0f),
36 register_int("matrix_size", &R);
37 register_float("circle_radius", &circle_radius);
38 register_float("gaussian_radius", &gaussian_radius);
39 register_float("correlation", &correlation);
40 register_float("noise", &noise);
43 string DeconvolutionSharpenEffect::output_fragment_shader()
46 sprintf(buf, "#define R %u\n", R);
49 assert(R <= 25); // Same limit as Refocus.
52 return buf + read_file("deconvolution_sharpen_effect.frag");
57 // Integral of sqrt(r² - x²) dx over x=0..a.
58 float circle_integral(float a, float r)
65 return 0.25f * M_PI * r * r;
67 return 0.5f * (a * sqrt(r*r - a*a) + r*r * asin(a / r));
70 // Yields the impulse response of a circular blur with radius r.
71 // We basically look at each element as a square centered around (x,y),
72 // and figure out how much of its area is covered by the circle.
73 float circle_impulse_response(int x, int y, float r)
76 // Degenerate case: radius = 0 yields the impulse response.
77 return (x == 0 && y == 0) ? 1.0f : 0.0f;
80 // Find the extents of this cell. Due to symmetry, we can cheat a bit
81 // and pretend we're always in the upper-right quadrant, except when
82 // we're right at an axis crossing (x = 0 or y = 0), in which case we
83 // simply use the evenness of the function; shrink the cell, make
84 // the calculation, and down below we'll normalize by the cell's area.
85 float min_x, max_x, min_y, max_y;
90 min_x = abs(x) - 0.5f;
91 max_x = abs(x) + 0.5f;
97 min_y = abs(y) - 0.5f;
98 max_y = abs(y) + 0.5f;
100 assert(min_x >= 0.0f && max_x >= 0.0f);
101 assert(min_y >= 0.0f && max_y >= 0.0f);
103 float cell_height = max_y - min_y;
104 float cell_width = max_x - min_x;
106 if (min_x * min_x + min_y * min_y > r * r) {
107 // Lower-left corner is outside the circle, so the entire cell is.
110 if (max_x * max_x + max_y * max_y < r * r) {
111 // Upper-right corner is inside the circle, so the entire cell is.
115 // OK, so now we know the cell is partially covered by the circle:
125 // The edge of the circle is defined by x² + y² = r²,
126 // or x = sqrt(r² - y²) (since x is nonnegative).
127 // Find out where the curve crosses our given y values.
128 float mid_x1 = (max_y >= r) ? min_x : sqrt(r * r - max_y * max_y);
129 float mid_x2 = sqrt(r * r - min_y * min_y);
130 if (mid_x1 < min_x) {
133 if (mid_x2 > max_x) {
136 assert(mid_x1 >= min_x);
137 assert(mid_x2 >= mid_x1);
138 assert(max_x >= mid_x2);
140 // The area marked A in the figure above.
141 float covered_area = cell_height * (mid_x1 - min_x);
143 // The area marked B in the figure above. Note that the integral gives the entire
144 // shaded space down to zero, so we need to subtract the rectangle that does not
145 // belong to our cell.
146 covered_area += circle_integral(mid_x2, r) - circle_integral(mid_x1, r);
147 covered_area -= min_y * (mid_x2 - mid_x1);
149 assert(covered_area <= cell_width * cell_height);
150 return covered_area / (cell_width * cell_height);
153 // Compute a ⊙ b. Note that we compute the “full” convolution,
154 // ie., our matrix will be big enough to hold every nonzero element of the result.
155 MatrixXf convolve(const MatrixXf &a, const MatrixXf &b)
157 MatrixXf result(a.rows() + b.rows() - 1, a.cols() + b.cols() - 1);
158 for (int yr = 0; yr < result.rows(); ++yr) {
159 for (int xr = 0; xr < result.cols(); ++xr) {
162 // Given that x_b = x_r - x_a, find the values of x_a where
163 // x_a is in [0, a_cols> and x_b is in [0, b_cols>. (y is similar.)
165 // The second demand gives:
167 // 0 <= x_r - x_a < b_cols
168 // 0 >= x_a - x_r > -b_cols
169 // x_r >= x_a > x_r - b_cols
170 int ya_min = yr - b.rows() + 1;
172 int xa_min = xr - b.rows() + 1;
175 // Now fit to the first demand.
176 ya_min = max<int>(ya_min, 0);
177 ya_max = min<int>(ya_max, a.rows() - 1);
178 xa_min = max<int>(xa_min, 0);
179 xa_max = min<int>(xa_max, a.cols() - 1);
181 assert(ya_max >= ya_min);
182 assert(xa_max >= xa_min);
184 for (int ya = ya_min; ya <= ya_max; ++ya) {
185 for (int xa = xa_min; xa <= xa_max; ++xa) {
186 sum += a(ya, xa) * b(yr - ya, xr - xa);
190 result(yr, xr) = sum;
196 // Similar to convolve(), but instead of assuming every element outside
197 // of b is zero, we make no such assumption and instead return only the
198 // elements where we know the right answer. (This is the only difference
200 // This is the same as conv2(a, b, 'valid') in Octave.
202 // a must be the larger matrix of the two.
203 MatrixXf central_convolve(const MatrixXf &a, const MatrixXf &b)
205 assert(a.rows() >= b.rows());
206 assert(a.cols() >= b.cols());
207 MatrixXf result(a.rows() - b.rows() + 1, a.cols() - b.cols() + 1);
208 for (int yr = b.rows() - 1; yr < result.rows() + b.rows() - 1; ++yr) {
209 for (int xr = b.cols() - 1; xr < result.cols() + b.cols() - 1; ++xr) {
212 // Given that x_b = x_r - x_a, find the values of x_a where
213 // x_a is in [0, a_cols> and x_b is in [0, b_cols>. (y is similar.)
215 // The second demand gives:
217 // 0 <= x_r - x_a < b_cols
218 // 0 >= x_a - x_r > -b_cols
219 // x_r >= x_a > x_r - b_cols
220 int ya_min = yr - b.rows() + 1;
222 int xa_min = xr - b.rows() + 1;
225 // Now fit to the first demand.
226 ya_min = max<int>(ya_min, 0);
227 ya_max = min<int>(ya_max, a.rows() - 1);
228 xa_min = max<int>(xa_min, 0);
229 xa_max = min<int>(xa_max, a.cols() - 1);
231 assert(ya_max >= ya_min);
232 assert(xa_max >= xa_min);
234 for (int ya = ya_min; ya <= ya_max; ++ya) {
235 for (int xa = xa_min; xa <= xa_max; ++xa) {
236 sum += a(ya, xa) * b(yr - ya, xr - xa);
240 result(yr - b.rows() + 1, xr - b.cols() + 1) = sum;
248 void DeconvolutionSharpenEffect::update_deconvolution_kernel()
250 // Figure out the impulse response for the circular part of the blur.
251 MatrixXf circ_h(2 * R + 1, 2 * R + 1);
252 for (int y = -R; y <= R; ++y) {
253 for (int x = -R; x <= R; ++x) {
254 circ_h(y + R, x + R) = circle_impulse_response(x, y, circle_radius);
258 // Same, for the Gaussian part of the blur. We make this a lot larger
259 // since we're going to convolve with it soon, and it has infinite support
260 // (see comments for central_convolve()).
261 MatrixXf gaussian_h(4 * R + 1, 4 * R + 1);
262 for (int y = -2 * R; y <= 2 * R; ++y) {
263 for (int x = -2 * R; x <= 2 * R; ++x) {
265 if (gaussian_radius < 1e-3) {
266 val = (x == 0 && y == 0) ? 1.0f : 0.0f;
268 val = exp(-(x*x + y*y) / (2.0 * gaussian_radius * gaussian_radius));
270 gaussian_h(y + 2 * R, x + 2 * R) = val;
274 // h, the (assumed) impulse response that we're trying to invert.
275 MatrixXf h = central_convolve(gaussian_h, circ_h);
276 assert(h.rows() == 2 * R + 1);
277 assert(h.cols() == 2 * R + 1);
279 // Normalize the impulse response.
281 for (int y = 0; y < 2 * R + 1; ++y) {
282 for (int x = 0; x < 2 * R + 1; ++x) {
286 for (int y = 0; y < 2 * R + 1; ++y) {
287 for (int x = 0; x < 2 * R + 1; ++x) {
292 // r_uu, the (estimated/assumed) autocorrelation of the input signal (u).
293 // The signal is modelled a standard autoregressive process with the
294 // given correlation coefficient.
296 // We have to take a bit of care with the size of this matrix.
297 // The pow() function naturally has an infinite support (except for the
298 // degenerate case of correlation=0), but we have to chop it off
299 // somewhere. Since we convolve it with a 4*R+1 large matrix below,
300 // we need to make it twice as big as that, so that we have enough
301 // data to make r_vv valid. (central_convolve() effectively enforces
302 // that we get at least the right size.)
303 MatrixXf r_uu(8 * R + 1, 8 * R + 1);
304 for (int y = -4 * R; y <= 4 * R; ++y) {
305 for (int x = -4 * R; x <= 4 * R; ++x) {
306 r_uu(x + 4 * R, y + 4 * R) = pow(double(correlation), hypot(x, y));
310 // Estimate r_vv, the autocorrelation of the output signal v.
311 // Since we know that v = h ⊙ u and both are symmetrical,
312 // convolution and correlation are the same, and
313 // r_vv = v ⊙ v = (h ⊙ u) ⊙ (h ⊙ u) = (h ⊙ h) ⊙ r_uu.
314 MatrixXf r_vv = central_convolve(r_uu, convolve(h, h));
315 assert(r_vv.rows() == 4 * R + 1);
316 assert(r_vv.cols() == 4 * R + 1);
318 // Similarly, r_uv = u ⊙ v = u ⊙ (h ⊙ u) = h ⊙ r_uu.
319 MatrixXf r_uu_center = r_uu.block(2 * R, 2 * R, 4 * R + 1, 4 * R + 1);
320 MatrixXf r_uv = central_convolve(r_uu_center, h);
321 assert(r_uv.rows() == 2 * R + 1);
322 assert(r_uv.cols() == 2 * R + 1);
324 // Add the noise term (we assume the noise is uncorrelated,
325 // so it only affects the central element).
326 r_vv(2 * R, 2 * R) += noise;
328 // Now solve the Wiener-Hopf equations to find the deconvolution kernel g.
329 // Most texts show this only for the simpler 1D case:
331 // [ r_vv(0) r_vv(1) r_vv(2) ... ] [ g(0) ] [ r_uv(0) ]
332 // [ r_vv(-1) r_vv(0) ... ] [ g(1) ] = [ r_uv(1) ]
333 // [ r_vv(-2) ... ] [ g(2) ] [ r_uv(2) ]
334 // [ ... ] [ g(3) ] [ r_uv(3) ]
336 // (Since r_vv is symmetrical, we can drop the minus signs.)
338 // Generally, row i of the matrix contains (dropping _vv for brevity):
340 // [ r(0-i) r(1-i) r(2-i) ... ]
342 // However, we have the 2D case. We flatten the vectors out to
343 // 1D quantities; this means we must think of the row number
344 // as a pair instead of as a scalar. Row (i,j) then contains:
346 // [ 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) ... ]
348 // g and r_uv are flattened in the same fashion.
350 // Note that even though this matrix is block Toeplitz, it is _not_ Toeplitz,
351 // and thus can not be inverted through the standard Levinson-Durbin method.
352 // There exists a block Levinson-Durbin method, which we may or may not
353 // want to use later. (Eigen's solvers are fast enough that for big matrices,
354 // the convolution operation and not the matrix solving is the bottleneck.)
356 // One thing we definitely want to use, though, is the symmetry properties.
357 // Since we know that g(i, j) = g(|i|, |j|), we can reduce the amount of
358 // unknowns to about 1/4th of the total size. The method is quite simple,
359 // as can be seen from the following toy equation system:
361 // A x0 + B x1 + C x2 = y0
362 // D x0 + E x1 + F x2 = y1
363 // G x0 + H x1 + I x2 = y2
365 // If we now know that e.g. x0=x1 and y0=y1, we can rewrite this to
367 // (A+B+D+E) x0 + (C+F) x2 = 2 y0
368 // (G+H) x0 + I x2 = y2
370 // This both increases accuracy and provides us with a very nice speed
372 MatrixXf M(MatrixXf::Zero((R + 1) * (R + 1), (R + 1) * (R + 1)));
373 MatrixXf r_uv_flattened(MatrixXf::Zero((R + 1) * (R + 1), 1));
374 for (int outer_i = 0; outer_i < 2 * R + 1; ++outer_i) {
375 int folded_outer_i = abs(outer_i - R);
376 for (int outer_j = 0; outer_j < 2 * R + 1; ++outer_j) {
377 int folded_outer_j = abs(outer_j - R);
378 int row = folded_outer_i * (R + 1) + folded_outer_j;
379 for (int inner_i = 0; inner_i < 2 * R + 1; ++inner_i) {
380 int folded_inner_i = abs(inner_i - R);
381 for (int inner_j = 0; inner_j < 2 * R + 1; ++inner_j) {
382 int folded_inner_j = abs(inner_j - R);
383 int col = folded_inner_i * (R + 1) + folded_inner_j;
384 M(row, col) += r_vv((inner_i - R) - (outer_i - R) + 2 * R,
385 (inner_j - R) - (outer_j - R) + 2 * R);
388 r_uv_flattened(row) += r_uv(outer_i, outer_j);
392 LLT<MatrixXf> llt(M);
393 MatrixXf g_flattened = llt.solve(r_uv_flattened);
394 assert(g_flattened.rows() == (R + 1) * (R + 1)),
395 assert(g_flattened.cols() == 1);
397 // Normalize and de-flatten the deconvolution matrix.
398 g = MatrixXf(R + 1, R + 1);
400 for (int i = 0; i < g_flattened.rows(); ++i) {
403 if (y == 0 && x == 0) {
404 sum += g_flattened(i);
405 } else if (y == 0 || x == 0) {
406 sum += 2.0f * g_flattened(i);
408 sum += 4.0f * g_flattened(i);
411 for (int i = 0; i < g_flattened.rows(); ++i) {
414 g(y, x) = g_flattened(i) / sum;
417 last_circle_radius = circle_radius;
418 last_gaussian_radius = gaussian_radius;
419 last_correlation = correlation;
423 void DeconvolutionSharpenEffect::set_gl_state(GLuint glsl_program_num, const string &prefix, unsigned *sampler_num)
425 Effect::set_gl_state(glsl_program_num, prefix, sampler_num);
429 if (fabs(circle_radius - last_circle_radius) > 1e-3 ||
430 fabs(gaussian_radius - last_gaussian_radius) > 1e-3 ||
431 fabs(correlation - last_correlation) > 1e-3 ||
432 fabs(noise - last_noise) > 1e-3) {
433 update_deconvolution_kernel();
435 // Now encode it as uniforms, and pass it on to the shader.
436 float samples[4 * (R + 1) * (R + 1)];
437 for (int y = 0; y <= R; ++y) {
438 for (int x = 0; x <= R; ++x) {
439 int i = y * (R + 1) + x;
440 samples[i * 4 + 0] = x / float(width);
441 samples[i * 4 + 1] = y / float(height);
442 samples[i * 4 + 2] = g(y, x);
443 samples[i * 4 + 3] = 0.0f;
447 set_uniform_vec4_array(glsl_program_num, prefix, "samples", samples, (R + 1) * (R + 1));