Implementing convolution from scratch

Specific Domains Signal and Image Processing
,
1

Hello,
I already posted this question on Stackoverflow, 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

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)

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, 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

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

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?

2

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

3

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.

4

Here is an implementation that works (except for borders)

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

5

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.

6

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:

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
1 Like