{"id":1495,"date":"2010-09-07T09:33:45","date_gmt":"2010-09-07T06:33:45","guid":{"rendered":"https:\/\/yehar.com\/blog\/?p=1495"},"modified":"2019-10-13T10:00:56","modified_gmt":"2019-10-13T07:00:56","slug":"phased-gaussian-convolution","status":"publish","type":"post","link":"https:\/\/yehar.com\/blog\/?p=1495","title":{"rendered":"Circularly symmetric convolution and lens blur"},"content":{"rendered":"

This article describes approaches for efficient isotropic two-dimensional convolution <\/strong>with disc-like and arbitrary circularly symmetric convolution kernels, and also discusses lens blur effects.<\/p>\n

Keywords: depth of field, circle of confusion<\/em>, bokeh<\/em>, circular blur<\/em>, lens blur<\/em>, hexagonal blur, octagonal blur<\/em>, real-time, DOF\n<\/em><\/p>\n

Gaussian function approach<\/h2>\n

The circularly symmetric 2-d Gaussian kernel is linearly separable<\/em>; the convolution can be split into a horizontal convolution followed by a vertical convolution. This makes 2-d Gaussian convolution relatively cheap, computationally, and \"Gaussian blur\" has become, partially for this reason, a popular operation in image processing. No other circularly symmetric (isotropic) convolution kernel is linearly separable. Usually a Gaussian kernel is real-valued. However, we can parameterize it with a complex number to obtain a complex 2-d Gaussian function that still satisfies conditions for separability and is circularly symmetric. This section of the article shows that the real part of weighted sums of such separable 2-d kernels can produce arbitrary 2-d circularly symmetric kernels. Special focus is given to implementation of circular blur, or convolution with a solid disk, which is useful for simulation of bokeh<\/em> (circle of confusion) of camera lenses with fully open apertures.<\/p>\n

\"\"<\/a>
A Gaussian convolution kernel, used in Gaussian blur (black = -maximum value, grey = 0, white = maximum value)<\/figcaption><\/figure><\/p>\n

A horizontal convolution followed by a vertical convolution (or in the opposite order) by 1-d kernels $f(x)$ and $f(y)$ effectively gives a 2-d convolution by 2-d kernel $f(x)\\times f(y)$. To produce separable circularly symmetric 2-d convolution, the condition that must be satisfied for all $(x, y)$ is $f(\\sqrt{x^2+y^2})\\times f(0) = f(x)\\times f(y)$. The right side of the equation is the 2-d kernel. The left side can be understood as taking a horizontal slice $f(x)\\times f(0)$ of the kernel at $y=0$ and substituting $x\\to r$ in the slice function to obtain a function $f(r)\\times f(0)$ of radius $r = \\sqrt{x^2+y^2}$ which the circularly symmetric kernel must be equal to at all $(x, y)$. If we assume that the 1-d kernels are normalized by $f(0)=1$, then the constraint for circular symmetry simplifies to $f(\\sqrt{x^2+y^2}) = f(x)\\times f(y)$. We will work with that simplifying assumption.<\/p>\n

For starters, let's treat magnitude and phase of the complex kernels separately. The circular symmetry condition may be broken down into a magnitude condition and a phase condition that must both be satisfied. The magnitude condition reads $|f(\\sqrt{x^2+y^2})| = |f(x) \\times f(y)| = |f(x)|\\times|f(y)| $, meaning that the kernel must have a Gaussian function as its magnitude function. Because in multiplication of complex numbers the phases are summed, the phase (argument) condition reads $\\theta(\\sqrt{x^2+y^2}) = \\theta(x)+\\theta(y)$, modulo $2\\pi$, with the shorthand $\\theta(x) = \\arg(f(x))$. The only family of phase functions to satisfy this is $\\theta(x) = a x^2$. The kernel function can be constructed as the product of a real Gaussian function, here called the envelope<\/em>, and a complex phasor<\/em> with $x^2$ as its argument. In one dimension these kernels have the form $e^{-a x^2} e^{i b x^2} = e^{- a x^2 + i b x^2} = e^{(- a + b i) x^2}$, where $i$ is the imaginary unit, $x$ is the spatial coordinate, $a$ is the spatial scaling constant of the envelope, and $b$ is the spatial scaling constant of the complex phasor. In two dimensions the form is $e^{- a x^2 + i b x^2} e^{- a y^2 + i b y^2} = e^{- a (x^2 + y^2) + i b (x^2 + y^2)} = e^{(- a + b i) (x^2 + y^2)}$, $y$ being the coordinate in the other dimension. These forms are not normalized to preserve average intensity in image filtering, but to give a value of 1 for $x = 0, y = 0$. Let's see what these look like.<\/p>\n\n\n\n\n
\"\"<\/a>
The real part of a phased Gaussian kernel<\/figcaption><\/figure><\/td>\n
\"\"<\/a>
The imaginary part of a phased Gaussian kernel<\/figcaption><\/figure><\/td>\n<\/tr>\n
\"\"<\/a>
The real part of a phased Gaussian kernel of infinite spatial scale of the magnitude part, showing the phasing only<\/figcaption><\/figure><\/td>\n
\"\"<\/a>
The imaginary part of the phased Gaussian kernel of infinite spatial scale of the magnitude part, showing the phasing only<\/figcaption><\/figure><\/td>\n
\"\"<\/a>
Magnitude of a phased Gaussian kernel, identical to the plain Gaussian kernel<\/figcaption><\/figure><\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n

Arbitrary circularly symmetric shapes can be constructed by a weighted sum of the real parts and imaginary parts of complex phasors of different frequencies. Fourier series is the standard approach for such designs. These give a repetitive pattern, but the repeats could be diminished by adjusting the spatial scale of the envelope. As necessary, sloping of the envelope can be compensated for in the weights of the sum, but simultaneous optimization of both will give the best results, as it can take advantage of that the envelopes need not be the same for each of the component kernels.<\/p>\n

Let's try a 5-component square wave approximation, designed manually with the Fourier series applet<\/a> of Paul Falstad.<\/p>\n\n\n\n
\"\"<\/a>
5-component equiripple square wave Fourier series<\/figcaption><\/figure><\/td>\n
\"\"<\/a>
Real part of a square wave approximated by 5 phased Gaussian components of infinite spatial scale of the magnitude part.<\/figcaption><\/figure><\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n

The formula for this series is $0.54100\\ +\\ 0.66318\\ \\cos(x)\\ -\\ 0.20942\\ \\cos(3 x)\\ +\\ 0.10471\\ \\cos(5 x)\\ -\\ 0.08726\\ \\cos(7 x)$. A bit of scaling (by 0.9) and shifting (by -0.01) was then done to make the ripples go evenly around the zero level and to avoid values greater than 1.<\/p>\n

This is a pretty good indicator of the capabilities of this approach. The weighted component kernels were:<\/p>\n\n\n\n
\"\"<\/a>
Component 0<\/figcaption><\/figure><\/td>\n
\"\"<\/a>
Component 1<\/figcaption><\/figure><\/td>\n
\"\"<\/a>
Component 2<\/figcaption><\/figure><\/td>\n
\"\"<\/a>
Component 3<\/figcaption><\/figure><\/td>\n
\"\"<\/a>
Component 4<\/figcaption><\/figure><\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n

By lowering the spatial scale of the magnitude part of each component kernel, one gets a result like this:<\/p>\n

\"\"<\/a>
First try on a circular convolution kernel<\/figcaption><\/figure><\/p>\n

Not so bad! The halos would need to be eliminated, and the disk is also fading a bit toward the edges.<\/p>\n

I searched for better circular kernels by global optimization<\/a> with an equiripple cost function, with the number of component kernels and the transition bandwidth as the parameters. Transition bandwidth defines how sharp the edge of the disk is. The sharper it is, the more ringing there will be both inside and outside the disk. A transition bandwidth setting of 1 means that the width of the transition is as large as the radius of the disk inside. Below, 1-d component kernel formulas are given for a pass band spanning $x \\in -1..1$ and a dual stop band spanning $ x\\ \\in\\ -1.2\\ ..\\ -1\\ \\cup\\ 1\\ ..\\ 1.2$.<\/p>\n\n\n\n
\"\"<\/a>
Number of components: 1, transition bandwidth: 0.200000, ripple: \u00b10.232417. Component 0: (cos(x*x*1.624835) * 0.767583 + sin(x*x*1.624835) * 1.862321) * exp(-0.862325*x* x)<\/figcaption><\/figure><\/td>\n
\"\"<\/a>
1-d view<\/figcaption><\/figure><\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n\n\n\n
\"\"<\/a>
Number of components: 2, transition bandwidth: 0.200000, ripple: \u00b10.075459. Component 0: (cos(x*x*5.268909) * 0.411259 + sin(x*x*5.268909) * -0.548794) * exp(-0.886528*x *x) Component 1: (cos(x*x*1.558213) * 0.513282 + sin(x*x*1.558213) * 4.561110) * exp(-1.960518*x* x)<\/figcaption><\/figure><\/td>\n
\"\"<\/a>
1-d view of the components and their sum (magenta)<\/figcaption><\/figure><\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n\n\n\n
\"\"<\/a>
Number of components: 3, transition bandwidth: 0.200000, ripple: \u00b10.026297. Component 0: (cos(x*x*5.043495) * 1.621035 + sin(x*x*5.043495) * -2.105439) * exp(-2.176490*x *x) Component 1: (cos(x*x*9.027613) * -0.280860 + sin(x*x*9.027613) * -0.162882) * exp(-1.019306* x*x) Component 2: (cos(x*x*1.597273) * -0.366471 + sin(x*x*1.597273) * 10.300301) * exp(-2.815110* x*x)<\/figcaption><\/figure><\/td>\n
\"\"<\/a>
1-d view of the components and their sum (cyan)<\/figcaption><\/figure><\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n\n\n\n
\"\"<\/a>
Number of components: 4, transition bandwidth: 0.200000, ripple: \u00b10.010843. Component 0: (cos(x*x*1.553635) * -5.767909 + sin(x*x*1.553635) * 46.164397) * exp(-4.338459* x*x) Component 1: (cos(x*x*4.693183) * 9.795391 + sin(x*x*4.693183) * -15.227561) * exp(-3.839993* x*x) Component 2: (cos(x*x*8.178137) * -3.048324 + sin(x*x*8.178137) * 0.302959) * exp(-2.791880*x *x) Component 3: (cos(x*x*12.328289) * 0.010001 + sin(x*x*12.328289) * 0.244650) * exp(-1.342190* x*x)<\/figcaption><\/figure><\/td>\n
\n

\"\"<\/a>
1-d view of the components and their sum (brown)<\/figcaption><\/figure><\/p>\n

\"\"<\/a>
Close-up<\/figcaption><\/figure><\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n\n\n\n
\"\"<\/a>
Number of components: 5, transition bandwidth: 0.200000, ripple: \u00b10.004062. Component 0: (cos(x*x*1.685979) * -22.356787 + sin(x*x*1.685979) * 85.912460) * exp(-4.892608 *x*x) Component 1: (cos(x*x*4.998496) * 35.918936 + sin(x*x*4.998496) * -28.875618) * exp(-4.711870 *x*x) Component 2: (cos(x*x*8.244168) * -13.212253 + sin(x*x*8.244168) * -1.578428) * exp(-4.052795 *x*x) Component 3: (cos(x*x*11.900859) * 0.507991 + sin(x*x*11.900859) * 1.816328) * exp(-2.929212* x*x) Component 4: (cos(x*x*16.116382) * 0.138051 + sin(x*x*16.116382) * -0.010000) * exp(-1.512961 *x*x)<\/figcaption><\/figure><\/td>\n
\n

\"\"<\/a>
1-d view of the components and their sum (magenta)<\/figcaption><\/figure><\/p>\n

\"\"<\/a>
Close-up<\/figcaption><\/figure><\/td>\n
\n

\"\"<\/a>
Close-up of the pass band<\/figcaption><\/figure><\/p>\n

\"\"<\/a>
Close-up of the stop band<\/figcaption><\/figure><\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n

The last one of the above, with 5 component kernels, is good enough for most practical applications as the ripple is as small as 1\/250, about as much as the least significant bit for 8-bits-per-channel graphics, so invisible in almost all situations. As the number of components is increased, the components have an increasingly large weighted amplitude, as much as 36 times their sum for the 5-component system. This may become a practical problem if one aims for convolution kernel shapes with sharp transitions. A construction via Fourier series would probably give better-behaving components, but require more components than abusing the Gaussian envelope for ripple control as seems to be the case in the above given composite kernels.<\/p>\n

Here's a bonus, with 6 components:<\/p>\n

\"\"<\/a>
Number of components: 6, transition bandwidth: 0.200000, ripple: 0.001918. Component 0: (cos(x*x*2.079813) * -82.326596 + sin(x*x*2.079813) * 111.231024) * exp(-5.14377 8*x*x) Component 1: (cos(x*x*6.153387) * 113.878661 + sin(x*x*6.153387) * 58.004879) * exp(-5.612426 *x*x) Component 2: (cos(x*x*9.802895) * 39.479083 + sin(x*x*9.802895) * -162.028887) * exp(-5.98292 1*x*x) Component 3: (cos(x*x*11.059237) * -71.286026 + sin(x*x*11.059237) * 95.027069) * exp(-6.5051 67*x*x) Component 4: (cos(x*x*14.810520) * 1.405746 + sin(x*x*14.810520) * -3.704914) * exp(-3.869579 *x*x) Component 5: (cos(x*x*19.032909) * -0.152784 + sin(x*x*19.032909) * -0.107988) * exp(-2.20190 4*x*x)<\/figcaption><\/figure><\/p>\n

The 1-d composite kernel can be truncated to range $x \\in -(\\textrm{pass band radius + transition bandwidth}\\ ..\\ \\textrm{pass band width + transition bandwidth})$. This only necessitates a stop band equiripple condition for $|x| \\in \\textrm{pass band radius + transition bandwidth}\\ ..\\ \\sqrt{2}(\\textrm{pass band radius+ transition bandwidth})$ to ensure a flat stop band in the corners of the 2-d convolution kernel square. However, at least for up to the above given 5th order composite kernel, the stop band ripples start to slope out already within that range, so no savings can be made.<\/p>\n

Let's see the 5-component circular blur in action.<\/p>\n\n\n\n\n\n\n
\"\"<\/a>
The standard test image, Lena<\/figcaption><\/figure><\/td>\n
<\/td>\n<\/tr>\n
\"\"<\/a>
Lena after 5-component circular blur with a disk diameter of 128 pixels and an additional transition band of 25.6 + 25.6 pixels . (contrast adjusted manually)<\/figcaption><\/figure><\/td>\n
\"\"<\/a>
Gaussian blur of similar size<\/figcaption><\/figure><\/td>\n<\/tr>\n
\"\"<\/a>
With halved blur disk radius, 32 pixels (contrast set manually)<\/figcaption><\/figure><\/td>\n
\"\"<\/a>
Gaussian blur of similar size<\/figcaption><\/figure><\/td>\n<\/tr>\n
\"\"<\/a>
With 16 pixels disk radius (contrast set manually)<\/figcaption><\/figure><\/td>\n
\"\"<\/a>
Gaussian blur of similar size<\/figcaption><\/figure><\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n

The circular blur here is the weighted sum of the real and imaginary parts of the 5 component convolutions.<\/p>\n

Let's try another image, one which will better show the blur disks.<\/p>\n\n\n\n\n
\"\"<\/a>
Hubble Deep View<\/figcaption><\/figure><\/td>\n
\"\"<\/a>
Hubble Deep View after 5-component circular blur (arbitrary contrast)<\/figcaption><\/figure><\/td>\n<\/tr>\n
\"\"<\/a>
With halved radius, 32 pixels<\/figcaption><\/figure><\/td>\n
\"\"<\/a>
disk radius halved once more, now 16 pixels (the bands at the top and at the bottom are due to some programming bug)<\/figcaption><\/figure><\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n

Looks like an out-of-focus camera, doesn't it?<\/p>\n

Recursive approximations<\/h3>\n

The examples presented above were done with 1-d finite impulse response<\/em> (FIR) filters of which I now have a FFT-based implementation. More efficient would be to use identical pairs of complex 1-d causal and anti-causal infinite impulse response<\/em> (IIR) filters to achieve a symmetric composite filter. Early results:<\/p>\n\n\n\n
\"\"<\/a>
First attempt at tail-truncated parallel poles approximation of phased Gaussian components of a 3-component composite kernel (on the right). 3 complex poles were used for each direction, employing the structure: (left+right)*(up+down) and identical coefficient sets for all directions. The coefficients used have not yet been re-optimized to increase the resulting kernel quality directly.<\/figcaption><\/figure><\/td>\n
\"\"<\/a>
The \"exact\" composite kernel. Note that the ringing outside the disc is suppressed in the approximation (on the left) due to tail truncation. The transition band spans -1.2..-1 and 1..1.2.<\/figcaption><\/figure><\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n

The above first attempt filter is quite heavy, about 110 complex multiplications (some of which involve real numbers) per color channel per pixel. But it is independent of the size of the blur (disregarding boundary conditioning that may add a penalty proportional to the radius times the number of edge pixels in the image).<\/p>\n

Perhaps circular blur is not the best use for the phased Gaussian kernels, as no sharp edges can be achieved. How about constructing an airy disc<\/em> for anisotropic low-pass filtering of the image in a manner that gives an equal upper frequency limit for features in diagonal, horizontal, vertical, or any direction? Food for thought... Update: yes, it works<\/a><\/p>\n

Updates<\/h3>\n

There are now GPU implementations of this method out there. Kleber Garcia gave a SIGGRAPH 2017 talk (accompanying conference paper<\/a>) about his 1 and 2-component blur implementations. They can be seen at work in FIFA 17 screenshots<\/a>. The talk inspired Bart Wronski to write a Shadertoy implementation<\/a>, with an accompanying text<\/a>. Kleber Garcia responded with his own Shadertoy implementation<\/a>.<\/p>\n

I have also gained more understanding of the math involved. Each 1-d component kernel is really just a Gaussian function, but with a complex width parameter. The original condition I wrote for separability demands that A) the 2-d convolution kernel is a circularly symmetrical revolution of the 1-d kernel (the value of the 2-d kernel at radius $r$ equals the 1-d kernel at position $r$), B) the same coefficients are used for the horizontal and the vertical 1-d kernels. Then the final output for each component is a weighted sum of the real and imaginary parts of the 2-d convolution (implemented as a horizontal and a vertical 1-d convolution) result. The weighed sum with a weight $a$ for the real part and a weight $b$ for the imaginary part can be seen as complex multiplication<\/a> by $a-bi$ followed by discarding of the imaginary part. If we modify requirement A to: \"the 2-d convolution kernel is circularly symmetrical\", then we can absorb the said multiplication into the 1-d convolution kernel by multiplying each coefficient by the complex number $\\sqrt{a-bi}$. Because we have two 1-d convolution passes, that number will be squared, hence the square root. Now we have coefficient-wise identical horizontal and vertical 1-d convolutions that produce a 2-d convolution with a kernel that is circularly symmetrical. The work saved is that we don't need to calculate (or store!) the imaginary output of the 2nd 1-d convolution, because only the real part of its result will be used.<\/p>\n

Having done this optimization, the 1-d kernel can be expressed (recycling variable names) simply as $(c+di)e^{(a + b i) x^2} = e^{ax^2}\\big(c\\cos(bx^2) - d\\sin(bx^2)\\big) + ie^{ax^2}\\big(d\\cos(bx^2) + c\\sin(bx^2)\\big),$ resulting in a 2-d kernel $(c+di)^2e^{(a+bi)(x^2+y^2)}.$ One can then globally optimize the coefficients $(a, b, c, d)$ of a horizontal slice $(c+di)^2e^{(a+bi)x^2}$ (note the square) of the 2-d kernel to get a desired shape of the real part $e^{ax^2}\\big((c^2 - d^2)\\cos(bx^2) - 2cd\\sin(bx^2)\\big)$ of the slice, or the sum of these from multiple components.<\/p>\n

Using this notation, the $(a, b, c, d)$ coefficients are obtained from the old coefficients $(a, b, r, i)$, where $r$ and $i$ are the weights of the real and the imaginary part of the component:<\/p>\n

a = -a;\r\nb = b;\r\nc = sqrt(2)*(sqrt(r*r + i*i) + r)*sqrt(sqrt(r*r + i*i) - r)\/(2*i);\r\nd = -sqrt(2)*sqrt(sqrt(r*r + i*i) - r)\/2;\r\n\r\n\/\/ Use these coefs as: \r\n\/\/ real_coef = (c*cos(x*x*a) - d*sin(x*x*a)) * exp(b*x*x);\r\n\/\/ imag_coef = (d*cos(x*x*a) + c*sin(x*x*a)) * exp(b*x*x);\r\n\r\n\/\/ a, b, c, d:\r\n\/\/ 1-component\r\n-0.8623250000, 1.6248350000, 1.1793828124, -0.7895320249\r\n\/\/ 2-component\r\n-0.8865280000, 5.2689090000, -0.7406246191, -0.3704940302\r\n-1.9605180000, 1.5582130000, 1.5973700402, -1.4276936105\r\n\/\/ 3-component\r\n-2.1764900000, 5.0434950000, -1.4625695191, -0.7197739911\r\n-1.0193060000, 9.0276130000, -0.1480093005, -0.5502424493\r\n-2.8151100000, 1.5972730000, 2.2293886172, -2.3101178772\r\n\/\/ 4-component\r\n-4.3384590000, 1.5536350000, 4.5141678065, -5.1132787901\r\n-3.8399930000, 4.6931830000, -3.7350649493, -2.0384600009\r\n-2.7918800000, 8.1781370000, 0.0866540887, -1.7480940853\r\n-1.3421900000, 12.3282890000, 0.3569701172, -0.3426757426\r\n\/\/ 5-component\r\n-4.8926080000, 1.6859790000, 5.7626795783, -7.4542110865\r\n-4.7118700000, 4.9984960000, -6.4033389291, -2.2547313456\r\n-4.0527950000, 8.2441680000, -0.2167382954, -3.6413223544\r\n-2.9292120000, 11.9008590000, 1.0940793322, -0.8300714338\r\n-1.5129610000, 16.1163820000, -0.3717954486, -0.0134482550\r\n\/\/ 6-component\r\n-5.1437780000, 2.0798130000, 5.2941931370, -10.5050024737\r\n-5.6124260000, 6.1533870000, 10.9927011254, -2.6383360349\r\n-5.9829210000, 9.8028950000, -10.1550051566, -7.9777845753\r\n-6.5051670000, 11.0592370000, 4.8737688428, -9.7488280697\r\n-3.8695790000, 14.8105200000, -1.6383505756, -1.1306841329\r\n-2.2019040000, 19.0329090000, -0.1309780866, -0.4122368969\r\n<\/pre>\n
Negative values in the 2-d convolution kernel may be off-putting, because they create a cheap-looking sharpening kind of effect. I optimized some low component count coefficient sets that give almost no negative values in the 2-d convolution kernel:<\/p>\n
\/\/ a, b, c, d\r\n\/\/ 1-component, stp-band maximum = 0.3079699786\r\n-2.0126256644, 1.0462519386, 1.7772360464, -1.5705215513\r\n\/\/ 2-component, stop-band maximum = 0.1422940279\r\n-1.8755934460, 5.6868281810, -1.1341154250, -0.4727562391\r\n-3.0859247607, 0.0311853116, 14.9843211153, -14.9911604837\r\n<\/pre>\n
\"\"<\/a>
Corner-to-corner slice of the (mostly) nonnegative 1 and 2-component 2-d kernels.<\/figcaption><\/figure><\/p>\n
As noted by Bart Wronski, in addition to solid discs, also donuts may be desirable bokeh shapes. Donuts are associated at least with large focal length mirror lenses. It would probably be possible to optimize the components to get a visually impressive sharp donut, sharper than the discs shown. I did some experiments about that, but it didn't look very nice.<\/p>\n
For approximation of general 2-d (or 3-d) kernels that do not need to be circularly symmetric, as sums of separable 2-d convolutions, see Treitel, S. & Shanks, J. (1971). The design of multistage separable planar filters.<\/a> IEEE Transactions on Pattern Analysis and Machine Intelligence, 9(1), 10\u201327.<\/p>\n
Circular blur using prefix sums and summed area tables<\/h2>\n
Circular blur done using the phased gaussian kernels is not perfect because it has the transition band. A completely different approach would be to calculate each blurred output pixel separately by summing the pixels that fall within the disk radius from the position of the pixel. There is a very efficient way to do this. The area of the disk can be subdivided into a minimum number of rectangular pieces and a summed area table<\/a> can be used to efficiently calculate the sum of pixels within each rectangle. Here is an illustration of the idea:<\/p>\n
\"\"<\/a>
A circle packed with a small number of rectangles.<\/figcaption><\/figure><\/p>\n
Kyle McDonald<\/a><\/strong> has been working on the same problem and shows a more efficient packing with his Processing application<\/a>:<\/p>\n
\"\"<\/a>
A more efficient packing by Kyle McDonald, requiring for this particular circle 60 summed area table lookups (or alternatively a total of 40 summed area table or sum table lookups).<\/figcaption><\/figure><\/p>\n
A somewhat simpler way would be to subdivide the circle into one-pixel-height strips. The sum of the pixel values of a strip could be calculated from the prefix sum over the length of each column, requiring only two lookups per strip. The number of rows would be roughly the diameter of the circle, and the number of lookups roughly four times the radius, or $4r$.<\/p>\n
By combining the use of summed area tables and horizontal and vertical prefix sums we get a rather boring-looking hybrid approach that appears to be the most efficient:<\/p>\n
\"\"<\/a>
A packing requiring less lookups than McDonald's if both summed area tables and row and column prefix sums are used. For this particular circle a total of 32 summed area table or prefix sum lookups are required.<\/figcaption><\/figure><\/p>\n
For large $r$, the hybrid approach takes about $(8-4\\sqrt{2})r \\approx 2.34\\ r$ lookups: 4 for the rectangle (from a summed area table) but dominated for large $r$ by 2 for each stripe (from horizontal or vertical prefix sum). In comparison, Kyle McDonald's approach takes about $(12-4\\sqrt{5})r \\approx 3.06\\ r$ table lookups with both summed area tables and prefix sums available and about $(16-\\frac{24\\sqrt{5}}{5})r \\approx 5.27\\ r$ table lookups with only the summed area table available. The latter result is identical for the hybrid approach. If the whole image is convolved with a disc of constant radius, the center rectangle in the hybrid approach can be done with a constant complexity box blur<\/a> and a summed area table will not be needed, relieving requirements for the bit depth used in the calculations.<\/p>\n
To get an anti-aliased edge for the circle, the edge pixels can be summed or subtracted with weights separately, doable in the ready-for-prefix-sum representation. The number of pixels that should receive anti-aliasing is about $(4\\sqrt{2} + 16 - \\frac{24\\sqrt{5}}{5})r \\approx 10.9\\ r$, increasing the total cost of the algorithm to $(24 - \\frac{24\\sqrt{5}}{5})r \\approx 13.27\\ r$. Each weight will be shared by 8 edge pixels, which may give some computational savings. In conclusion, anti-aliased or not, the complexity of this way of blurring will be proportional to the radius of the disc times the number of pixels in the image. The decompositions (lookup recipes) of different size blurs could be stored in tables, which would allow also the use of a depth map<\/em> to blur different parts of the image with a different disk size.<\/p>\n
Alternatively to exact antialising, one could approximate the longer monotonic antialias gradations as linear, by use of two horizontal prefix sum passes in series and two vertical prefix sum passes in series. The slope of each linear section would be thrown in first, then after calculation of the first prefix sum, the constant term of the linear section would be added, along with a value that terminates the section. Either of these would coincide with the values that define where the flat section of the disc starts for that horizontal or vertical stripe. The method could be made more accurate by increasing the number of linear sections, perhaps up to two per monotonous gradation.<\/p>\n
If different size blurs are used in the same image, then one decision must be made on how the result should look: whether the blur scatters<\/em> or gathers<\/em>. Let's demonstrate with a square box blur, increasing the blur size towards the top of the image:<\/p>\n\n\n\n
\"\"<\/a>
Input image for box blur, just a single bright pixel on black background. Size of the box blur will vary by row as indicated by the width of the blue area<\/figcaption><\/figure><\/td>\n
\"\"<\/a>
Scattering box blur with a width of 5, as designated for the position of the bright input pixel, scatters the intensity of the pixel to the surrounding 5 times 5 pixel area.<\/figcaption><\/figure><\/td>\n
\"\"<\/a>
Gathering box blur resulting from outputting at each position the average intensity of surrounding input pixels from an area the size of which is dictated by the output position.<\/figcaption><\/figure><\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n
With a variable blur kernel, scattering blur is more correct, because it preserves average image intensity, unlike gathering blur. With smooth blur size gradations, gathering blur should also look okay.<\/p>\n
If the kernel is constant, blurring can be described mathematically by convolution<\/a> of the input image by the blur kernel: $o = i*K$, where $o$ is the output image, $i$ is the input image and $K$ is the blur kernel. There are two ways to directly implement the convolution: In the scattering way, first the output image is cleared to zero, and then for each input pixel, a copy of the kernel $K$ is multiplied by the intensity of the input pixel, spatially shifted to the position of the input pixel, and added to the output image. In the gathering way each output pixel is formed by shifting a copy of the kernel reversed along each axis, to the position of the output pixel, and setting the output pixel to the sum of products of the pixel intensities in the shifted and reversed kernel and the input image. If the same kernel is used throughout the image, both ways of doing convolution will be identical. But if the size of the blur varies, we must acknowledge whether we do a scattering or a gathering convolution by $K$. We denote these by $\\hat{*}$ for a scattering and $\\check{*}$ for a gathering convolution by the variable kernel. The mnemonic is that the arrowhead points away to indicate scattering.<\/p>\n
We define the order of calculation of written formulas to be from left to right with the exception that multiplication and convolution are calculated first and addition and subtraction last. The convolution operator $*$ has the following algebraic properties<\/a> :<\/p>\n