How do I create an SVector/SMatrix with a function and arbitrary indices?

Question

The prof demonstrated a 2D Gaussian kernel using an SMatrix (see pic below). I am trying to implement a 1D kernel of arbitrary size in a similar fashion. I have tried code as below and am unable to find any good resource online on the correct way to initialize it. All the examples are static with very static elements at the point of initializing.

begin
    G_1D(i, σ) = ℯ^(-(i^2)/(2*(σ^2)))/sqrt(2*π*(σ^2))
    G_2D(x, y, σ) = ℯ^(-((x^2)+(y^2))/(2*(σ^2)))/sqrt(2*π*(σ^2))
end

begin
    function gaussian_kernel_1D(n, σ = 1)
        if iseven(n)
            mid = n÷2
            kernel = SVector{n}(map(i -> G_1D(i, σ)), CartesianIndex(-mid:mid-1))
            kernel ./ sum(kernel)
            return kernel
        else
            mid = n÷2
            kernel = SVector{n}(map(i -> G_1D(i, σ)), CartesianIndex(-mid:mid))
            kernel ./ sum(kernel)
            return kernel
        end
    end
    
    function gaussian_kernel_2D(n, σ = 1)
        if iseven(n)
            mid = n÷2
            kernel = SMatrix{n,n}(map(x,y -> G_2D(x,y, σ)), CartesianIndices((-mid:mid-1, -mid:mid-1)))
            kernel ./ sum(kernel)
            return kernel
        else
            mid = n÷2
            kernel = SMatrix{n,n}(map(x,y -> G_2D(x,y, σ)), CartesianIndices((-mid:mid-1, -mid:mid-1)))
            kernel ./ sum(kernel)
            return kernel
        end
    end
end

1. Is using an SVector the right solution? If so, could you help me with the right way to initialize it?
2. If the SVector is not the right idea, what would you recommend?

Here's the error detail on calling the gaussian kernel function.

MethodError: no method matching (::Main.workspace68.var"#2#6"{Int64})()
Closest candidates are:
#2(!Matched::Any) at /Users/
map(::Main.workspace68.var"#2#6"{Int64})@abstractarray.jl:2247
gaussian_kernel_1D(::Int64, ::Int64)@Other: 10
gaussian_kernel_1D(::Int64)

I'm a relatively new learner, so do point out if the code is incorrect, it will be very helpful.

Answer 1

Your problem isn't with SVector, but with the way you call G_1D, more specifically, you get an error here:

julia> CartesianIndex(-mid:mid)
ERROR: MethodError: no method matching CartesianIndex(::UnitRange{Int64})

You could use CartesianIndices, but I don't see why you would. Just use -mid:mid directly, with broadcasting. like this:

G_1D.(-mid:mid, σ)

So then

kernel = SVector{n}(G_1D.(-mid:mid, σ))

For 2D you can do this:

kernel = SMatrix{n,n}(G_2D.(-mid:mid, (-mid:mid)', σ))

There's another problem, though, here:

kernel ./ sum(kernel)
return kernel

This does not change the value of kernel. You are basically doing the same as this:

kernel_nonnormalized = kernel
kernel_normalized = kernel ./ sum(kernel)
return kernel_nonnormalized

If you just write

return kernel ./ sum(kernel)

then you get a normalized value out.

Whether you are justified in using StaticArrays is a bit hard to say. The dimensions should be small, and ideally their sizes should be known at compile time, while in your code it is determined at runtime. That may not be a problem though, as long as this happens at an isolated place, for example at the beginning of you calculation. Just try it out and see.

Edit: BTW, if you don't mind my saying so, you are violating the DRY (Don't Repeat Yourself) principle unnecessarily, and missing out on using multiple dispatch. Here's one way to avoid that (may not be the optimal one):

gauss(i, σ) = ℯ^(-(i^2)/(2*(σ^2)))/sqrt(2*π*(σ^2))
gauss(x, y, σ) = ℯ^(-((x^2)+(y^2))/(2*(σ^2)))/sqrt(2*π*(σ^2))

function makekernel(f, n, σ=1)
    ind = (0:n-1) .- (n÷2)
    kernel = SVector{n}(f.(ind, σ))
    return kernel ./ sum(kernel)
end

function makekernel(f, n::NTuple{2}, σ=1)
    ind1 = (0:n[1]-1) .- (n[1]÷2)
    ind2 = (0:n[2]-1) .- (n[2]÷2)
    kernel = SMatrix{n[1], n[2]}(f.(ind1, ind2', σ))
    return kernel ./ sum(kernel)
end

You just need one name for the gauss function, no need for separate names for 1D and 2D versions. Also makekernel can have the same name in both cases. No need to create separate function names with 1D and 2D inside the name.

Notice also that you don't need separate code paths for even and odd lengths.

Here's how you would call them:

kernel1D = makekernel(gauss, 7)  # here's 1D kernel
kernel2D = makekernel(gauss, (7, 6))  # here's a 7x6 2D kernel

Now there's also room for using different functions than gauss, and you also see a pattern for extending this to higher dimensions.