The box blur gave quite a bit of artifacts. Even a larger kernel will introduce artifacts. The main problem is that the filter is strong in the x- and y-directions.

A way to avoid this is to use a filter that is radially symmetric. A common such filter with that property is a Gaussian filter. We multiply the input image with a 2D Gaussian function that depends on the distance.

A Gaussian function usually has three free parameters:

  • $\sigma$, which denotes the size of the kernel
  • $A$, which denotes the amplitude of the kernel
  • $b$, which denotes the center

We want the center of the kernel to be on the current pixel. The parameter $\sigma$ will remain a free parameter, which we can adjust freely:

In [1]:
import halide as hl

import imageio
from halide_book.interactive import Slider, interact
from halide_book.paths import Paths
import halide.imageio

sigma = hl.Param(hl.Float(32), "sigma", 2.0)
sigma_slider = Slider(sigma, min_value=0.1, max_value=4.0, step=0.01)

The usual way to make the kernel normalized is to set $A = (\sigma\sqrt{2\pi})^{-1}$. However, we might not be able to make it infinitely large. Thus, we should rather normalize based on the sum of the values. We therefore use $A = 1$ in the definition of the kernel, and proceed to calculate the sum of all it components, which we subsequently divide by to obtain a normalized kernel.

The Gaussian kernel is defined as

In [2]:
x = hl.Var("x")
y = hl.Var("y")
c = hl.Var("c")

gaussian = hl.Func("gaussian")
gaussian[x, y] = hl.exp(-0.5 * (x**2 + y**2) / sigma**2)
gaussian.compute_root()
Out[2]:
<halide.Func 'gaussian'>

We can plot the kernel on its own and see how sigma affects the size of the Gaussian:

In [3]:
interact(gaussian, [(-8, 17), (-8, 17)], controls=[sigma_slider])

sigma: 2.0

The normalization factor then becomes

In [4]:
rdom = hl.RDom([(-3, 7), (-3, 7)])
normalization_factor = hl.Func("normalization_factor")
normalization_factor[()] = 1.0 / hl.sum(gaussian[rdom.x, rdom.y])

normalization_factor.compute_root()
Out[4]:
<halide.Func 'normalization_factor'>

Next, we define our blurred output image as a convolution

In [5]:
image_data = halide.imageio.imread(Paths.static / "stuff.png")
image_buffer = hl.Buffer(image_data)
bounded_buffer = hl.BoundaryConditions.constant_exterior(image_buffer, exterior=0)

input_values = hl.Func("input_values")
input_values[x, y, c] = hl.f32(bounded_buffer[x, y, c] / 255.0)

gaussian_blurred = hl.Func("gaussian_blurred")
gaussian_blurred[x, y, c] = 0.0
gaussian_blurred[x, y, c] += (
    normalization_factor[()]
    * input_values[
        x + rdom.x,
        y + rdom.y,
        c,
    ]
    * gaussian[
        rdom.x,
        rdom.y,
    ]
)

blurred = hl.Func("blurred")
blurred[x, y, c] = hl.u8(gaussian_blurred[x, y, c] * 255.0)

input_values.compute_root()
gaussian_blurred.compute_at(blurred, x)

interact(
    blurred,
    controls=[sigma_slider],
    size=[800, 600, 3],
)

sigma: 2.0