# Implementing convolution from scratch

**URL:** <https://discourse.julialang.org/t/implementing-convolution-from-scratch/93723>\
**Category:** Signal and Image Processing\
**Tags:** image-processing, convolution\
**Created:** [January 29, 2023, 7:03am UTC](https://discourse.julialang.org/t/implementing-convolution-from-scratch/93723 "2023-01-29T07:03:32Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![skittishdev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/skittishdev/32/46340_2.png) [@skittishdev](https://discourse.julialang.org/u/skittishdev)\
**Post date:** [January 29, 2023, 7:03am UTC](https://discourse.julialang.org/t/implementing-convolution-from-scratch/93723/1 "2023-01-29T07:03:32Z")

</div>

Hello,  
I already posted this question on [Stackoverflow](https://stackoverflow.com/questions/75271900/implementing-convolution-from-scratch-in-julia), but decided I’d ask here as well.

Please find my question from Stackoverflow pasted below (verbatim).

I am trying to implement convolution by hand in Julia. I’m not too familiar with image processing or Julia, so maybe I’m biting more than I can chew.

Anyway, when I apply this method with a 3\*3 edge filter `edge = [0 -1 0; -1 4 -1; 0 -1 0]` as `convolve(img, edge)`, I am getting an error saying that my values are exceeding the allowed values for the RGBA type.

### Code

```julia
function convolve(img::Matrix{<:Any}, kernel)
	(half_kernel_w, half_kernel_h) = size(kernel) .÷ 2
	(width, height) = size(kernel)
	cpy_im = copy(img)
	for row ∈ 1+half_kernel_h:height-half_kernel_h
		for col ∈ 1+half_kernel_w:width-half_kernel_w
			from_row, to_row = row .+ (-half_kernel_h, half_kernel_h)
			from_col, to_col = col .+ (-half_kernel_h, half_kernel_h)
			cpy_im[row, col] = sum((kernel .* RGB.(img[from_row:to_row, from_col:to_col])))
		end
	end
	cpy_im
end

```

### Error (original)

```julia-auto
ArgumentError: element type FixedPointNumbers.N0f8 is an 8-bit type representing 256 values from 0.0 to 1.0, but the values (-0.0039215684f0, -0.007843137f0, -0.007843137f0, 1.0f0) do not lie within this range.
See the READMEs for FixedPointNumbers and ColorTypes for more information.

```

I am able to identify a simple case where such error may occur (a white pixel surrounded by all black pixels or vice-versa). I tried “fixing” this by attempting to follow the advice [here from another stackoverflow question](https://stackoverflow.com/questions/53027923/how-do-you-handle-negative-pixel-values-after-filtering), but I get more errors to the effect of `Math on colors is deliberately undefined in ColorTypes, but see the ColorVectorSpace package.`.

### Code attempting to apply solution from the other SO question

```julia-auto
function convolve(img::Matrix{<:Any}, kernel)
	(half_kernel_w, half_kernel_h) = size(kernel) .÷ 2
	(width, height) = size(kernel)
	cpy_im = copy(img)
	for row ∈ 1+half_kernel_h:height-half_kernel_h
		for col ∈ 1+half_kernel_w:width-half_kernel_w
			from_row, to_row = row .+ [-half_kernel_h, half_kernel_h]
			from_col, to_col = col .+ [-half_kernel_h, half_kernel_h]
			cpy_im[row, col] = sum((kernel .* RGB.(img[from_row:to_row, from_col:to_col] ./ 2 .+ 128)))
		end
	end
	cpy_im
end

```

### Corresponding error

```julia-auto
MethodError: no method matching +(::ColorTypes.RGBA{Float32}, ::Int64)
Math on colors is deliberately undefined in ColorTypes, but see the ColorVectorSpace package.

Closest candidates are:
+(::Any, ::Any, !Matched::Any, !Matched::Any...) at operators.jl:591
+(!Matched::T, ::T) where T<:Union{Int128, Int16, Int32, Int64, Int8, UInt128, UInt16, UInt32, UInt64, UInt8} at int.jl:87
+(!Matched::ChainRulesCore.AbstractThunk, ::Any) at ~/.julia/packages/ChainRulesCore/a4mIA/src/tangent_arithmetic.jl:122

```

Now, I can try using `convert` etc., but when I look at the big picture, I start to wonder what the idiomatic way of solving this problem in Julia is. And that is my question. If you had to implement convolution by hand from scratch, what would be a good way to do so?

---

<div class="post-metadata">

**Author:** ![skittishdev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/skittishdev/32/46340_2.png) [@skittishdev](https://discourse.julialang.org/u/skittishdev)\
**Post date:** [January 29, 2023, 8:09am UTC](https://discourse.julialang.org/t/implementing-convolution-from-scratch/93723/2 "2023-01-29T08:09:49Z")

</div>

> [@skittishdev](#):
>
> `(width, height) = size(kernel)`

Found a small bug in this line. This has to be `(width, height) = size(img)`

---

<div class="post-metadata">

**Author:** ![GunnarFarneback](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gunnarfarneback/32/1827_2.png) [@GunnarFarneback](https://discourse.julialang.org/u/GunnarFarneback)\
**Post date:** [January 29, 2023, 8:21am UTC](https://discourse.julialang.org/t/implementing-convolution-from-scratch/93723/3 "2023-01-29T08:21:40Z")

</div>

`N0f8` is not the easiest type to do calculations with. I would recommend that you look for a way to convert your image to floating point before trying to convolve it. I would also suggest starting with gray scale and getting that to work before involving colors.

---

<div class="post-metadata">

**Author:** ![skittishdev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/skittishdev/32/46340_2.png) [@skittishdev](https://discourse.julialang.org/u/skittishdev)\
**Post date:** [January 29, 2023, 9:24am UTC](https://discourse.julialang.org/t/implementing-convolution-from-scratch/93723/4 "2023-01-29T09:24:43Z")

</div>

Here is an implementation that works (except for borders)

```julia
function convolve(img::Matrix{<:Any}, kernel)
	(half_kernel_h, half_kernel_w) = size(kernel) .÷ 2
	(height, width) = size(img)
	cpy_im = copy(img)
	# println(Dict("width" => width, "height" => height, "half_kernel_w" => half_kernel_w, "half_kernel_h" => half_kernel_h, "row range" => 1+half_kernel_h:(height-half_kernel_h), "col range" => 1+half_kernel_w:(width-half_kernel_w)))
	for row ∈ 1+half_kernel_h:(height-half_kernel_h)
		for col ∈ 1+half_kernel_w:(width-half_kernel_w)
			from_row, to_row = row .+ (-half_kernel_h, half_kernel_h)
			from_col, to_col = col .+ (-half_kernel_w, half_kernel_w)
			vals = Dict()
			for method ∈ [red, green, blue, alpha]
				x = sum((kernel .* method.(img[from_row:to_row, from_col:to_col])))
				if x > 1
					x = 1
				elseif x < 0
					x = 0
				end
				vals[method] = x
			end
			cpy_im[row, col] = RGBA(vals[red], vals[green], vals[blue], vals[alpha])
		end
	end
	cpy_im
end

```

This implementation does appear to have much poorer performance than `imfilter` though

---

<div class="post-metadata">

**Author:** ![skittishdev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/skittishdev/32/46340_2.png) [@skittishdev](https://discourse.julialang.org/u/skittishdev)\
**Post date:** [January 29, 2023, 9:29am UTC](https://discourse.julialang.org/t/implementing-convolution-from-scratch/93723/5 "2023-01-29T09:29:33Z")

</div>

Thank you @GunnarFarneback. Your comment made me want to split processing of RGBA separately, and I was able to come up with an implementation that works, though it is not the most performant.

---

<div class="post-metadata">

**Author:** ![GunnarFarneback](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gunnarfarneback/32/1827_2.png) [@GunnarFarneback](https://discourse.julialang.org/u/GunnarFarneback)\
**Post date:** [January 29, 2023, 9:43am UTC](https://discourse.julialang.org/t/implementing-convolution-from-scratch/93723/6 "2023-01-29T09:43:23Z")

</div>

There are many tricks for making these kinds of algorithms fast. This is an okay and mostly straightforward baseline for a zero-padded imfilter in arbitrary dimensions:

```julia
using OffsetArrays
function imfilter_zero_pad!(out::Array{T, N},
                            img::Array{T, N},
                            kernel::Array{T, N}) where {T, N}
    kernel = OffsetArrays.centered(kernel)
    kernel_indices = CartesianIndices(kernel)

    Threads.@threads for I in CartesianIndices(img)
        @inbounds begin
            s = zero(T)

            @simd for dI in kernel_indices
                J = I + dI
                if J ∈ CartesianIndices(img)
                    s += img[J] * kernel[dI]
                end
            end

            out[I] = s
        end
    end

    return
end

```