Re-posted from: http://learningjulia.com/2017/02/24/blurring-and-manipulation.html

Now that I’ve been able to load an image and play with some colors, I wanted to implement some basic image manipulation functionality: blurring. To handle out-of-bounds errors when back-projecting to the input image, I learned about image overloading, played with the Interpolations.jl package, and generated some really trippy mirror art.

Here is my notebook, annotated with my findings and pretty pictures.

If you want, you can skip to some headers below:

- Generating kernels
- Basic blurring
- Interpolations.jl
- Mirror interpolation
- Aside: broadcasting linear interpolations

In this notebook I want to play with image manipulation I’m going to implement a blurring function on an input image.

One fundamental concept in computational photography: iterate over the output to generate the appropriate input to use. The reason for this is that when you iterate over the output, you can compute which pixels from the input image(s) are used to compute the output. This ensures that every pixel in the output will be assigned a value.

So, to blur an image:

- Generate a kernel
- Allocate a new image (blurring can’t be an in-place operation because an input pixel is used multiple times, so can’t be overwritten)
- Iterate over the output, compute the input section used for that output, and convolve each input section with the kernel.

I am also curious to see the effects of different methods on the timing of the entire program.

```
using Images, TestImages, Colors, FileIO;
```

```
img = testimage("lighthouse")
```

## Generating kernels¶

I want to generate a few custom types (like structs in C++ or enum types in Python 3) to make data manipulation a little bit easier.

I tried to make a type called `size`

, but got this error message (I took a screenshot because I didn’t actually want to overwrite a new method called size). It looks like there is another method defined in the images packages – actually MANY of them! No matter, I will define a new type with a different name!

```
load("03-overwrite-main.jpg")
```

```
length(methods(size))
```

```
type coord
x::Int
y::Int
end
```

```
methods(coord)
```

An interesting feature of Julia types is that they define default function constructors for a type, hence the output of the `methods`

function above.

In the process of creating a function to generate a kernel in 2 dimensions, I want to generate a 1D gaussian. I can do this with an anonymous function. Even cooler, variable names in Julia can be represented as unicode! As a friend says, “it can lead you to do some really awesome things, and some things you probably should not do.” I think using LaTeX symbols in the Julia REPL is quite magical. Certainly easier to represent math!

(To generate unicode as a variable name, type `\`

followed by the appropriate LaTeX symbol and press tab.)

```
# interestingly, pi is an internal constant!
# and you can even use LaTeX to represent it
println(pi)
println(π)
println(e)
```

The function name is `\phi`

. This is the Gaussian function centered at 0, and you specify an input value `x`

and a standard deviation `\sigma`

.

```
ϕ(x, σ) = e^(-x^2 / (2 * σ^2)) / (σ * sqrt(2 * π))
```

The fact that I can use the Julia dot notation and have it broadcast the function appropriately is magical!

```
# generate a 1D gaussian kernel
ϕ.(-3:3, 4)
```

My coordinate type is initialized as an Int, but I want it to do something even better: I want to, when initializing a coordinate, if it is initialized as a float, have the coordinates rounded and turned into an `Int`

. I can take advantage of the fact that types define constructors by adding a default constructor inside the type definition.

```
type coord
x::Int
y::Int
coord(x, y) = new(Int(round(x)), Int(round(y)))
end
```

```
methods(coord)
```

```
# you can see that it works!
coord(2.0, 2.0)
```

```
function generate_gaussian_kernel(σ=3.0, kernel_size=coord(3, 3))
# because of the magic of functions and types, this just works!
center = coord(kernel_size.x / 2, kernel_size.y / 2);
kernel = Array{Float64}(kernel_size.x, kernel_size.y);
for x in 1:kernel_size.x
for y in 1:kernel_size.y
distance_to_center = sqrt((center.x - x)^2 + (center.y - y)^2);
kernel[x, y] = ϕ(distance_to_center, σ);
end
end
# normalize the coefficients at the end
# so the sum is 1
kernel /= sum(kernel);
# and assert it
@assert sum(kernel) - 1 < 1e-5;
# explicitly return a kernel
return kernel;
end
```

```
generate_gaussian_kernel(1, coord(3, 3))
```

## Putting it all together: blurring¶

Now let’s put it all together to generate a blurred version of the image using our kernels!

The first thing you need to do is make sure your input image is floats, so you can do proper arithmetic. After much Googling I found that you need to use the `convert`

function from Images.jl. Alternatively, you can use casting.

```
[convert(Array{RGB{Float64}}, img) RGB{Float64}.(img)]
```

```
function blur(img; kernel_size=15)
# explicitly allocate an array of zeroes, to avoid
# reading random memory, and ensuring the output
# image starts out black
img_blurred = zeros(RGB{Float64}, size(img, 1), size(img, 2));
kernel = generate_gaussian_kernel(1, coord(kernel_size, kernel_size));
for x in 1:size(img_blurred, 1) - kernel_size
for y in 1:size(img_blurred, 2) - kernel_size
# convolve the kernel with the image segment.
# the .* operator is element-wise multiplication
# as opposed to the * operator that is
# matrix multiplication
img_blurred[x, y] = sum(RGB{Float64}.(
@view img[x:x+kernel_size-1, y:y+kernel_size-1]) .* kernel)
end
end
return img_blurred;
end
```

```
img_blurred = blur(img)
```

## Interpolating nearest neighbors with the Interpolations package¶

The black pixels on the edges of the image don’t have exactly kernel-size matches in the input image, so for now I am just not not computing those pixels. The way I do this is by modifying the `for`

loop to not iterate over the last pixels. The right way to fix this is to implement an “out of bounds” function to access pixels outside the range of an image.

In fact, the Julia Interpolations.jl package already implemented this!

```
# Pkg.add("Interpolations");
```

```
using Interpolations;
```

You can set an extrapolation object on an interpolation object.

```
?extrapolate
```

```
# here I am using Constant interpolation and Constant
# (for some reason called Flat) extrapolation
# a.k.a. Nearest Neighbor interpolation
img_interp = extrapolate(interpolate(img, BSpline(Constant()), OnGrid()), Flat());
```

```
function blur_nn(img; kernel_size=15)
img_blurred = zeros(RGB{Float64}, size(img, 1), size(img, 2));
kernel = generate_gaussian_kernel(1, coord(kernel_size, kernel_size));
for x in 1:size(img_blurred, 1)
for y in 1:size(img_blurred, 2)
# convolve the kernel with the image segment.
# the .* operator is element-wise multiplication
# as opposed to the * operator that is
# matrix multiplication
img_blurred[x, y] = sum(RGB{Float64}.(
img_interp[x:x+kernel_size-1, y:y+kernel_size-1]) .* kernel)
end
end
return img_blurred;
end
```

```
img_blurred = blur_nn(img)
```

Interestingly enough, the last few rows in the blurred image are black. The reason for this that if you look at the last rows in the original image, it is a fully-black row:

```
img[end, :]
```