如何创建具有函数和任意索引的SVector / SMatrix？

问题描述

教授演示了使用SMatrix的2D高斯核（请参见下图）。我正在尝试以类似的方式实现任意大小的一维内核。我尝试了以下代码，但无法以正确的方式在线上找到任何好的资源对其进行初始化。所有示例都是静态的，在初始化时具有非常静态的元素。

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,CartesianIndices((-mid:mid-1,-mid:mid-1)))
            kernel ./ sum(kernel)
            return kernel
        else
            mid = n÷2
            kernel = SMatrix{n,-mid:mid-1)))
            kernel ./ sum(kernel)
            return kernel
        end
    end
end

1. 使用SVector是否是正确的解决方案？如果是这样，您能帮我用正确的方式初始化它吗？
2. 如果SVector不是正确的主意，您会推荐什么？

这是有关调用高斯内核函数的错误详细信息。

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)

我是一个相对较新的学习者，因此请指出如果代码不正确，将非常有帮助。

解决方法

您的问题不是SVector，而是您调用G_1D的方式，更具体地说，您在这里遇到错误：

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

您可以使用CartesianIndices，但我不知道为什么要这么做。只需将-mid:mid直接用于广播即可。像这样：

G_1D.(-mid:mid,σ)

那么

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

对于2D，您可以执行以下操作：

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

不过，这里还有另一个问题：

kernel ./ sum(kernel)
return kernel

这不会更改kernel的值。您基本上是这样做的：

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

如果你只是写

return kernel ./ sum(kernel)

然后您得到归一化的值。

您是否有理由使用StaticArrays很难说。尺寸应较小，理想情况下，尺寸应在编译时已知，而在代码中，尺寸应在运行时确定。但是，只要这发生在孤立的地方（例如在计算开始时），这可能就不是问题。只需试用一下即可。

编辑：顺便说一句，如果您不介意我的话，那就是在不必要地违反了DRY（请勿重复自己）原则，并且错过了使用多个分派的机会。这是避免这种情况的一种方法（可能不是最佳方法）：

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

对于高斯函数，您只需要一个名称，对于1D和2D版本，则不需要单独的名称。在这两种情况下，makekernel也可以具有相同的名称。无需在名称内使用1D和2D创建单独的函数名称。

还要注意，对于偶数和奇数长度，您不需要单独的代码路径。

这是您的称呼方式：

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

现在还可以使用与gauss不同的功能，并且您还会看到一种将其扩展到更高尺寸的模式。