Asymmetric recursive Gaussian filtering for space-

variant artificial bokeh

Tuan Q. Pham
Canon Information Systems Research Australia
Building A, The Park Estate, 5 Talavera Road,
Macquarie Park, NSW, 2113, Australia
tuan.pham@cisra.canon.com.au

Abstract—This paper describes a recursive filter that space-variant IIR filters in artificial bokeh application. Section
handles discontinuous space-variant blur. The two-dimensional V summarises the main contributions of this works.
space-variant blur is approximated by separable Gaussian
filtering along the x- and y-dimensions. Within each dimension,
II. RECURSIVE GAUSSIAN FILTERS
the one-dimensional Gaussian filter is further approximated by
a pair of forward and backward Infinite Impulse Response (IIR)
filters applied in parallel. By modifying the blur sigma’s of the A. Infinite impulse response filters
IIR filters as they approach a blur discontinuity, smearing There are two types of filters used in digital signal
artefacts are eliminated. Different modifications of the forward processing: FIR and IIR filters. FIR filters are implemented as
and backward filter sigma’s results in an asymmetric space- a weighted average of the input samples. The set of weights or
variant IIR Gaussian filter. Our asymmetric IIR filters produce filter coefficients forms the impulse response of the filter. Due
visually pleasing background blur for a complicated scene to a finite number of filter coefficients, this type of filters has
whose contents are at different depths to the camera. a finite impulse response. The FIR filters are often
characterised by its filter coefficients. For examples, [0.25 0.5
Keywords—artificial bokeh, IIR filter, recursive filter, space- 0.25] is a 3-tap low-pass filter, and [-1 0 1] is a 3-tap central
variant filter
gradient filter.

I. INTRODUCTION Unlike FIR filters, IIR filters have a feedback component

from a previously filtered output. Figure 1a illustrates this
Image filtering is a frequently used operation in digital concept in a discrete time domain, where the filter output y[n]
image processing. It is a building block of many applications at a time instance n is a linear combination of the input sample
including edge detection, image anti-aliasing, image resizing, x[n] and a feedback component y[n-1] from a previous time
etc. This paper presents an image enhancement application of instance n-1:
filtering called artificial bokeh. In artificial bokeh, a variable
amount of blur is added to each pixel in the image to create a = + . −1 (1)
gradual blur of the background. The foreground subjects,
which stay sharp, become the main focus of the image. This IIR filter, which exhibits an exponential decay
behaviour for 0<α<1, can also be characterised by its transfer
Digital signal filtering is usually referred to as a function in the Z domain [6]:
convolution of an input signal with a filter to produce an
output signal. The size and shape of the filter directly
( ) 1 (2)
correspond to the amount of blur applied to the input signal at ( )= =
any given point. This type of filtering has a Finite Impulse ( ) 1−
Response (FIR), which is a sampled version of the filter itself. The z-1 block in Figure 1b is a unit delay, which
The computational complexity of FIR filters is therefore corresponds to the delayed feedback component α.y[n-1] in
proportional to the size of the filter. While it is suitable for the discrete time domain. Due to the feedback, the impulse
small filters, FIR filtering is often too costly for large filters response of an IIR filter is infinite, that is given an impulse
that it has to be speeded-up in Fourier domain. input, an infinite number of non-zero values will come out (see
However, there is a more efficient way of image filtering Figure 1a). Because the impulse responses of IIR filters are
whose complexity is independent of the filter size. Infinite not truncated, IIR filters have much better frequency response
Impulse Response (IIR) filters achieve this using a fixed than FIR filters given the same number of coefficients.
number of filtering taps by feeding previously filtered output Another advantage of IIR filters is that through recursion they
pixels to the filtering at a current pixel. Up to now, this type use fewer taps to obtain an equivalent FIR filter of any filter
of recursive filters is suitable for smooth and slowly varying size.
blur only. We propose to extend existing IIR filters to handle
sudden variation in the blur map.
The remainder of this paper is organised as follows.
Section II reviews the theory of IIR filtering and different
approximations of the Gaussian filter using IIR filters. A
straightforward extension of existing IIR filters to handle a) exponential decay filter: 0<α<1 b) its block diagram in Z domain
space-variant blur is presented in Section IIC together with its
limitations. Our proposed modifications to the existing IIR Figure 1: A recursive filter in discrete time domain and Z domain.
filters to overcome their limitations are presented Section III.
Section IV demonstrates the suitability of our asymmetric

978-1-5386-6602-9/18/$31.00 ©2018 IEEE

 IIR filters are not without drawbacks, though. IIR filters where
can be unstable. For example, the exponential filter in Figure
1 shoots to infinity for | α |≥1. The infinite impulse response =
= ( − 1) (8)
can also result in unwanted signal bleeding over an extended
duration. = ( + 1)
=− (9)

B. Gaussian filters from IIR filters = −2 =2 (10)

The IIR filters presented in Section IIA are causal filters,
The filters from equation 5 and 6 are called second order
that is the filter output at time instance n depends only on the
recursive filters because the feedback delay is up to two
input samples at time instances before n. The impulse
samples before (or after) the current sample (n±2 versus n).
responses of such filters are zeros prior to the impulse. The
impulse response of a Gaussian filter, on the other hand, is 2) Gaussian filter from fourth-order IIR filters
symmetric with non-zeros components on both sides of the The Deriche’s second-order recursive filter in equation 3
central impulse. To achieve this symmetric impulse response, does not approximate the Gaussian function very well because
a combination of causal and anti-causal filters is often of a small number of filter coefficients. By increasing the filter
employed. Anti-causal filters are implemented similar to order to four, Farnebäck and Westin obtained a much better
causal filters except that the input samples are processed in approximation [5]. They used the following trigonometric
reverse order, starting from the last sample instance to the first exponential function to approximate a Gaussian function with
sample instance. σ=1 pixel:
Two ways to combine causal and anti-causal filters to
produce a symmetric filter response are depicted in Figure 2. ℎ( ) = ∑ ( cos + sin ) , t≥0 (11)
The filters can be applied sequentially as in Figure 2a or they The fit parameters ai, bi, ωi, and λi (i=1,2) was obtained by
can be applied in parallel followed by a summation as in Nelder-Mead simplex optimisation (fminsearch function
Figure 2b. In this paper, we focus on the parallel applications in Matlab).
of causal and anti-causal filters because it is easier to
manipulate the shape of the combined filter’s impulse The causal and anti-causal filters are fourth-order
response. recursive filters:

= + + + (12)
− − − −
= −∞, … , ∞

= + + + (13)
a) cascaded filter [12] b) parallel filter [4] − − − −
= ∞, … , −∞
Figure 2: Two schemes of combining causal and anti-causal filters
= + (14)
where the filter coefficients nj, mj, and dj (j=1...4) are
1) Gaussian filter from second-order IIR filters derived from the function fit parameters ai, bi, ωi, and λi (i=1,2)
Deriche is the first to propose using two independent and the scale σ. Details of this derivation can be found in the
passes of IIR filters to approximate the Gaussian filter [4]. The original paper [5] and the author’s implementation in ITK1.
impulse response of the approximated filter is:
Figure 3 compares Deriche’s second-order IIR filter and
( ) = ( | | + 1) | | (3) Farnebäck’s fourth-order IIR filter with the Gaussian filter at
scale σ=10. The impulse response of the fourth-order filter
where = 5⁄(2√ ) and (red dotted line) closely matches the Gaussian function (green
thick line) within four standard deviation range (-4σ ≤ t ≤ 4σ).
(1 − ) (4) Root Mean Squared Error (RMSE) of the fourth-order filter
= against the Gaussian function is two orders of magnitude
1+2 −
smaller than the RMSE of the second-order filter.
This filter can be realised from a summation of a causal
and anti-causal IIR filter:

= + − − (5)
∀ = −∞, … , ∞
= + − − (6)
∀ = ∞, … , −∞
= + (7)

1
Function itkRecursiveGaussianImageFilter.txx in the Insight
Segmentation and Registration Toolkit (ITK) is available at
http://lmi.bwh.harvard.edu/~gunnar/itk_recursive_gaussian.html
 -1
both causal and anti-causal directions at each pixel. While this
assumption is acceptable for slowly varying blur, it is not
0.05 10

suitable around blur discontinuities because the blur widths

0.04 -2
10

are different on either side of the pixel. Experimental results

0.03 -3
10
h(t)0.02 h(t)
10
-4 in Figure 4 confirm this, where the foveation blur in Figure 4a
Gaussian Gaussian σ=10
looks reasonable, but artefacts appear at abrupt blur transitions
0.01
nd
2 order 10
-5
2nd order, RMSE=2×10-3
0
4 nd
order -6
4nd order, RMSE=4×10-5 in Figure 4b. The artefacts are due to colour bleeding from the
-0.01
-40 -20 0
t
20 40
10
-40 -20 0
t
20 40 sharper region into the blurrier region.
a) linear response scale b) log response scale

Figure 3: Impulse responses of second and fourth order IIR

Gaussian filters.

3) Separable implementation for optimal cache usage

The Gaussian filters presented so far are one-dimensional
(1D) filters. When processing multi-dimensional (nD) images,
the 1D filters can be applied separably along each dimension
because the Gaussian function is separable. A common way a) foveation blur: 0≤σ≤20 b) discontinuous blur σ∈{0,10}
to do this is to first process the input image in the row
direction, followed by processing the output of the first step in Figure 4: Outputs of fourth-order IIR Gaussian filters for different
the column direction. Because the image is often stored in space-variant sigmas (the corresponding σ is shown in the inset).
memory and cache in a row-major order, row processing can
be very efficient. Processing along the column direction, Colour bleeding across blur discontinuities by a two-pass
however, involves continual read/write access to multiple (causal + anti-causal) IIR filters can be explained with the help
cache lines. This can result in frequent cache thrashing for of Figure 5. In this diagram, the input signal x is a step edge
large images, which slows down the processing. To prevent with collocated transition with the blur function σ. Each
cache thrashing, the image is transposed after it has been piecewise flat region in the input corresponds to a different
filtered in the first dimension [9]. Filtering the next dimension sigma: σ=0 when x=1, and σ=1 when x=0. The output is
will therefore happen along the rows of the transposed image. therefore expected to be unchanged. During a backward pass,
This improves data cache usage for both write-through and where the blur transits from one to zero, the anti-causal output
write-back caches. −
is indeed the same as the input x (Figure 5b). However,
during a forward pass, where the blur transits from zero to one,
C. Space-variant recursive Gaussian filter the causal filter (in red) with σ=1 for pixels next to the edge
In space-variant Gaussian filtering, the smoothing scale σ integrates over a region of mixed intensities on both side of
varies over the image. By treating the IIR filters’ coefficients the edge (orange region). This cause an intensity leakage from
as variables of space, Tan et al. used recursive filters for the sharp region (σ=0) into the blurry region (σ=1), and the
foveation blur [10]. However, this space-variant Gaussian shape of + changes accordingly. This intensity leakage can
smoothing only works for slowly varying σ. This is because be minimised by our proposed solution in the next section.
the feedback samples y[n-i] (i=1,2,…) should have been
blurred at roughly the same scale as the current output sample
y[n].
Space-variant recursive Gaussian filtering has also been
used for edge-preserving image denoising by Alvarez et al.
[1]. Motivated by the Perona-Malik anisotropy diffusion
model [8], Alvarez et al. proposed to vary the parameter α in
a) causal pass: bleeding b) anti-causal pass: no bleeding
the second-order IIR Gaussian filter (equation 3) according to
the gradient magnitude of the input signal x:
Figure 5: Intensity bleeding around a sudden blur discontinuity.

15
=ℎ ( )

where h(s) = Ms + ε is a positive function of s with M>0 III. ASYMMETRIC SPACE-VARIANT RECURSIVE GAUSSIAN
and ε=0.001. Because σ is inversely proportional to α, the FILTER
space-variant blur is small close to edges and large in flat We propose an asymmetric space-variant recursive filter
image areas. To improve stability, Alvarez et al.’s solution to stop intensity leakage at blur discontinuities. Instead of
requires three cascaded filters for each pass over the signal having the same blur σ for both the causal and anti-causal pass,
(causal and anti-causal), two of which are second-order IIR filters in each direction get a different σ. For 2D images, where
filters. For 2D images, Alvarez et al. averaged two separable recursive filtering is applied in two directions (forward and
filter results: x followed by y and y followed by x to minimize backward) and separably in two dimensions (x and y), all four
the influence of filtering order. sigmas , , , can be different. As the sliding filter
approaches a blur discontinuity along its path, the blur width
1) Limitations of existing space-variant IIR filters
can be reduced to minimise intensity leakage from the other
Both Tan et al.’s [10] and Alvarez et al.’s filters [1]
assume locally isotropic blur because the same σ is used for
side of the discontinuity (this applies when σ changes abruptly
from low to high as in Figure 5a).
0
A. Constraining the rate of change of the directional
sigmas
0

The idea of constraining the rate of change of the

directional sigmas to minimise intensity leakage across blur
discontinuities is illustrated in Figure 6. For the causal
direction, intensity leakage occurs when the blur sigma + 0
changes abruptly from low to high. To minimise this leakage,
the rate of increase of + is tapered to at most 1/3 per pixel.
The value of 1/3 was chosen based on the 3-sigma rule of the
Gaussian function which states that 99.7% of the area under
the function lie within three standard deviation. As illustrated 0
in the fourth row of Figure 6, the blur width of the causal filter
0

a) symmetric filter b) asymmetric filter

just after a rising σ is reduced so that the filter support does
not extend over the other side of the blur discontinuity.
Similarly, intensity leakage occurs along the anti-causal Figure 7: Impulse responses of a symmetric and an asymmetric
forward + backward IIR filters.
direction when the blur sigma − changes abruptly from high
to low. To minimise this leakage, the rate of decrease of − is
also tapered to at most 1/3 per pixel. The weighted identity filter can be further decomposed
into two parts wδ = + δ + − δ, where + = -n0(σ+)/2 and −
The different treatments of the directional sigmas around = n0( − )/2 (n0 is the value of the first recursive tap in the
blur discontinuities lead to asymmetric Gaussian filters. These forward direction in equation 5 or 12). n0(σ)/2 represents half
asymmetric filters avoid smoothing across blur the height of the filter response at the central tap. As a result,
discontinuities. subtracting n0(σ)δ/2 from the forward filter response and
1 adding n0(σ)δ/2 to the backward filter response make these
σ 0 0 one-sided filter responses symmetric about the central tap.
+ 1 Each of the filter responses after the delta correction sums up
σ 0 0
to 0.5.
̶ 1
σ 0 0 The complete signal flow diagram for asymmetric space-
causal filters anti-causal filters variant recursive Gaussian filtering is shown in Figure 8. After
1 1 the delta correction step, the forward and backward filter
x σ+=0 σ+=.1 σ ̶ =1 σ+= 1 σ̶ =.5 σ+=0 results at each pixel are combined using a weighted sum,
σ̶ =0 0 σ̶ =0 where the weights are proportioned to the corresponding blur
1 Σ Σ 1 sigma. This weighted combination was found to further
y 0 reduce halo artefacts around blur discontinuities. The halo was
present, albeit negligible, because IIR filters have infinitely
Figure 6: Asymmetric Gaussian filters avoid smoothing across long response that extends to the other side of the blur
blur discontinuities. discontinuity.

B. Asymmetric filter construction

The IIR filtering scheme in Section IIB does not produce
mirrored causal and anti-causal filters. The causal filter has a
non-zero response at the central tap, whereas the anti-causal
filter does not. This is due to the different input delays in the
causal and anti-causal IIR filter in equation 5-6 and 12-13. The
sum of these one-sided filter responses, however, form a
normalised symmetric Gaussian filter (see Figure 7a).
Unfortunately, the sum of a causal and anti-causal filter
response at different sigmas does not result in a normalised Figure 8: Asymmetric IIR filters stop bleeding across blur
filter. As a result, we have to add a weighted identity filter wδ discontinuities.
to the mix to make the total response of the asymmetric filter
equal to one (see Figure 7b).
Note that similar to existing space-variant IIR filters in
Section IIC, the filter coefficients ni, mi, and di are space-
variant. To avoid recomputing these recursive filter
coefficients at each pixel in real-time, their values can be pre-
computed, stored and retrieved from a look-up table [10].
Figure 9 demonstrates the results of our asymmetric
recursive Gaussian filter on the same image and space-variant
blur in Figure 4b. Note how the directional sigma’s in Figure
9a (where red=1, blue=0, and yellow=0.5) has been modified The regions pointed to by the arrows in Figure 11 are
from the input blur σ to reduce the rate of change at the rising enlarged in Figure 10 for side-by-side comparison. For each
and falling step edges. Unlike the result of Tan et al.’s filter region of interest, eight images are shown: the input image,
[10] in Figure 4b, our filter result in Figure 9b show no visible Deriche’s 2-tap filter result, Farnebäck’s 4-tap filter result,
artefacts at blur discontinuities. colour-coded blur map (top row), mip-map interpolation, our
2-tap and 4-tap asymmetric IIR filter results, and a filtered grid
+
σx
-
σx output. Mip-map interpolation is another space-variant
Gaussian blur method [11] computed by interpolation across
the scale axis of a Gaussian pyramid [3]. The 10-pixel-pitch
grid point image (superimposed on the colour coded blur map)
+ - was filtered using our 2-tap IIR filter to show the effective
σy σy
blurring kernel at each grid point.
Figure 10a show how our filters can render sharp
foreground hair against a blurry background without depth
boundary artefacts (when compared to the symmetric IIR
a) four directional sigmas b) discontinuous blur σ∈{0,10} filters) or colour bleeding (when compared to mip-map
interpolation). A closer look at the grid filtered image at the
Figure 9: Outputs of our asymmetric recursive Gaussian filter for bottom right confirms this: the foreground grid points stay
the same space-variant blur in Figure 4b (the blur σ is shown in the sharp, whereas the background grid points are completely
inset). smeared out given its blur σ equal the sampling pitch of the
grid.
IV. APPLICATIONS Figure 10b-c show regions with fine structures at different
depths from the camera, all of which undergo some blurring.
A. Artificial bokeh Once again, colour bleeding from sharper regions into
adjacent blurrier regions is evident in the mip-map
Images captured by SLR (single-lens reflex) cameras have
interpolated image. The 4-tap filters seem to produce more
a desirable characteristic of a narrow depth of field. The
artefacts than the 2-tap filters. This applies to both prior art’s
subject of focus is in maximum sharpness whereas the
symmetric filters and our asymmetric filters (Figure 10b). The
background, usually at a different depth to the camera, appears
grid filtered image show isotropic blur everywhere except at
blurred. The visually pleasing background blur, also known as
blur boundaries, where the effective blur kernel respects depth
bokeh, is only achievable optically under large apertures. As
edges.
a result, images from low-cost compact cameras often appear
sharp everywhere. This extended depth of field effect may
render the image unprofessional.
To make images from point-and-shoot cameras appear
more professional, artificial bokeh can be added to the original
sharp image. IIR Gaussian filter is perfectly suitable for this
purpose because it can efficiently apply variable amount of
blur to different areas in the image. To achieve a realistic
bokeh, one needs to have depth information at every pixel and
convert this depth to blur sigma. The depth information can be
obtained from a co-located depth sensor or it can be estimated
from the captured images using computational techniques
such as depth from defocus [7]. a) hair region (blue arrow in Figure 11, σblue=0, σred=10)
We illustrate the effectiveness of our asymmetric IIR
Gaussian filter in bokeh rendering given a full-resolution
depth map. Figure 11a shows a 1.5-megapixel input image
with false colours corresponding to an annotated segmentation
map. A single blur sigma is assigned to each of the six
coloured segments. The colours are coded so that cooler
colours corresponding to depths closer to the camera and
warmer colours for depths further away from the camera.
After assigning blur sigmas ranging from 0 to 10 to each
of the six different depths, artificial bokeh were synthesized in
Figure 11b. Four IIR filters were tested: the space-variant
versions of the 2-tap and 4-tap IIR filters from Deriche [4] and
b) table corner (red arrow in Figure 10, σgreen=4, σpink=8, σred=10)
Farnebäck [5] respectively and our corresponding asymmetric
filters in Section III. Both prior art filters produce severe
artefacts around depth discontinuities on the left of Figure 11b
where the red and blue arrows point. Our filters, on the other
hand, show little artefact while achieving pleasing background
blur.
 this type of artefacts around thin hair strands hovering in front
of a blurred background. Layer compositing [11] is a standard
approach to this problem because background layers at
different depths are blurred independently before being
blended together. For the case in Figure 13, background
behind the hair strands will be smoothly blurred across while
ignoring the thin foreground pixels. Note that our space-
variant blur filters can still be applied to each layer to simulate
the subtle change in blur across a face as often happens in a
portrait photo.

B. Edge-stopping blur
c) flowers (green arrow in Figure 10, σcyan=2, σyellow=6, σred=10)
A prior art bokeh rendering method uses anisotropic
Input Deriche‘s 2- Farnebäck’s 4- Blur map and point diffusion to perform space-variant edge-stopping blur [2].
tap IIR filter tap IIR filter grid Diffusion-based approaches produce good edge-stopping
N/A 11.4 24.2 N/A bokeh but they have one disadvantage: slow speed due to their
iterative nature. The proposed asymmetric recursive Gaussian
Mip-map Our 2-tap Our 4-tap
interpolation asymmetric IIR asymmetric IIR
Our 2-tap filtering of filter can be used as a large-step speedup of anisotropic
the point grid
filtering filtering diffusion. The diffusivity parameter in anisotropic diffusion is
0.487 12.6 25.2
5.28 c=g(‖∇I‖), where g is a function that maps [0 +∞) gradient
magnitude ‖∇I‖ values to [1 0] output values such as
d) Image order and run-time on Matlab 2012 for a 2MP image (in seconds) g ∇ = exp (− ∇ / ). This diffusivity parameter can be
used to directly control the blur width of the space-variant
Gaussian filter. In the function mapping, k is a pre-determined
Figure 10: Zoom-in on the regions of interests in Figure 11. constant, for example, k=median(‖∇I‖) is the median of the
gradient magnitude of the image. An example of the space-
Figure 10d tabulates the run-times of different space- variant diffusivity c for the Lena image is given in the inset of
variant filters under investigation. The 2-tap filters are two- Figure 14b. The diffusivity is small (dark intensity) around
time faster than the 4-tap filters due fewer number of filtering edges in the scene, and large (bright intensity) in
taps. Our asymmetric IIR filters are slightly slower than prior homogeneous image areas. By setting the blur width to a linear
art’s symmetric IIR filters because four different directional multiplication of this diffusivity, e.g. σ=16×c, the result of the
sigmas have to be computed from the input blur sigma. asymmetric recursive Gaussian filter closely resembles the
Filtering on a grayscale grid-point image is also more than filtered result of anisotropic diffusion (see Figure 14). The
two-time faster than on a 3-channel RGB image. Note that the gradient measure ‖∇I‖ can be computed from the input
mip-map interpolation method is significantly faster than the image itself or it can be computed from a guide image.
IIR filters. This is because our implementation of the non- Typically, the guide image is in good alignment with the input
vectorisable IIR filters is in pure Matlab code. We expect a image to allow preserving the same edges in both images. The
significantly faster run-time of the IIR filters if implemented guide image could be in a different modality compared to the
in C. To illustrate this point, we used a C/MEX input image. For example, the input image and the guide
implementation of space-invariant IIR filter in Matlab image could be a non-flash/flash image pair, or a
(filter.m) to filter the same 2MP colour image with a fixed visible/infrared image pair.
blur sigma of 5.7 pixel (average of the space-variant sigma
used in Figure 11) and measured a run-time of 0.349 seconds.
This is even faster than our mip-map implementation. Because
the complexity of IIR filters is independent of the blur sigma,
we expect this run-time figure is extendable to our asymmetric
space-variant IIR filters.
Similar to most post-processing bokeh synthesis methods,
our asymmetric recursive filters require a good depth map for
a visually pleasing background blur. The depth discontinuities
should ideally coincide with edges in the scene. However, a
small misalignment of the depth image is acceptable as seen a) 20 iterations of Perona-Malik b) our edge-stopping blur with
in Figure 12d. Artefacts only become visible when the anisotropic diffusion with λ=0.25,
(‖∇ ‖) = exp(−‖∇ ‖/5) = 16 × exp(−‖∇ ‖)
misalignment of depth is significant compared to the size of
object of interest as in Figure 12f. This type of artefacts due to
depth map error can be alleviated using error concealment Figure 14: Anisotropic diffusion versus space-variant Gaussian blur.
[11].
A key differentiation of our space-variant recursive blur One potential drawback of approximating anisotropic
method over previous methods is also a cause of undesirable diffusion by a single iteration of separable Gaussian filtering
artefacts in certain bokeh synthesis cases. Because our filters in the X-direction followed by a single iteration of separable
completely block colour diffusion across depth Gaussian filtering in the Y-direction is directional dependent
discontinuities, background blur may not appears natural on artefacts for large blur width. Different orders of filtering (X
either side of a thin foreground structure. Figure 13 illustrates followed by Y or Y followed by X) result in slightly different
outputs. The discrepancy due to filtering order is more REFERENCES
pronounced for larger blur width. In order to reduce these [1] L. Alvarez, R. Deriche, and F. Santana, “Recursivity and PDE’s in
directional dependent artefacts, the desired blur can be image processing”, in Proc. of ICPR, pp. 242-248, 2000.
implemented with multiple passes of X- followed by Y- [2] M.Bertalmio, R.Fort, and D.Canchez-Crespo, “Real-time, accurate
filtering. Usually, 2 iterations of separable filtering are enough depth of field using anisotropic diffusion and programmable graphics
to reduce these artefacts to a negligible level. A 2D Gaussian cards”, in Proc. of 3DPVT, 2004.
filter with blur width σ can be implemented with two [3] P.J. Burt and E.H. Adelson, “The Laplacian pyramid as a compact
applications of a separable Gaussian filter with blur width image code”, IEEE Trans. on Communication 31(4): 532-540, 1983.
⁄√2. [4] R. Deriche, “Fast algorithms for low-level vision”, PAMI, 12(1): 78-
87, 1990.
[5] G. Farnebäck and C.-F. Westin, “Improving Deriche-style recursive
V. CONCLUSION AND FUTURE WORKS Gaussian filters”, Journal of Mathematical Imaging and Vision,
26:293-299, 2006.
We have presented a new algorithm for space-variant
[6] A.V. Oppenheim and R.W. Schafer, Discrete-Time Signal Processing,
Gaussian blur using edge-stopping recursive filters. The new 2nd Edition, Prentice Hall, 1999.
algorithm constrains the blur spread around blur [7] A.P. Pentland, “A new sense for depth of field”, PAMI, 9(4):523-531,
discontinuities to avoid colour bleeding across the 1987.
discontinuity. This is especially useful in artificial bokeh [8] P. Perona and J. Malik, “Scale-space and edge detection using
application where blur sigma can change abruptly across anisotropic diffusion”, PAMI, 12(7):629-639, 1990.
depth boundaries. [9] K. Schutte and G.M.P. van Kempen, “Optimal cache usage for
separable image processing algorithms on general purpose
Compared to existing IIR Gaussian filters, our asymmetric workstations”, Signal Processing, 59(1):113-122, 1997.
space-variant IIR Gaussian filters produce significantly less [10] S. Tan, J.L. Dale, and A. Johnston, “Performance of three recursive
artefacts at the cost of being slightly slower. algorithms for fast space-variant Gaussian filtering”, Real-Time
Imaging, 9:215-228, 2003.
For future works, we plan to re-implement our filtering [11] S. Webster, D. Katchalov, and T.Q. Pham, “Variable blend width
algorithm in C using OpenMP to fully utilise the parallel compositing”, US patent 9311901, filed on 27 June 2014.
processing cores in modern CPUs. We expect the C [12] I.T. Young and L.J. van Vliet, “Recursive implementation of the
implementation will deliver sub-second run-time for colour Gaussian filter”, Signal Processing, 44:139-151, 1995.
images at HDTV resolution (i.e. 2 mega-pixel RGB image).

a) 1442×1080 image with 2 depth planes (σred=0, σblue=10) b) synthesised background blur image using our 2-tap
filter

c) <4-pixel depth map mis- d) acceptable bokeh given e) depth map misalignment f) noticeable artefacts
alignment around the wig small depth boundary error around flower petals (~10 pixels) due to large depth error

Figure 12: Bokeh synthesis artefacts due to a misaligned depth map.

a) 1442×1080 input image with 2 depth planes b) synthesised background blur image c) thin hair artefacts

Figure 13: Bokeh synthesis artefacts around thin foreground structures (2-pixel-thick hair strands).

a) input image with depth planes in false colour (σblue=0, σcyan=2, σgreen=4, σyellow=6, σpink=8 (spiky plant behind the
guitar), σred=10)

Deriche's 2-tap symmetric IIR Our 2-tap asymmetric IIR

Farneback's 4-tap symmetric IIR Our 4-tap asymmetric IIR

b) Artificial bokeh rendering using four different space-variant IIR filters

Figure 11: Bokeh synthesis using a high-quality depth map.