Q: FiniteDifferences | gradient of a 3D scalar array

Specific Domains Numerics
1

Hello,

I have a 3D array consisting of numerical scalar values:

imgVlm::Array{Float64, 3}

The grid spacing along each axis is uniform, but the spacing values are different for each axis.

FiniteDifferences seems to be the package to use as opposed to the basic

diff(A::AbstractArray; dims::Integer)

From the FiniteDifference, somewhat sparse, documentation, my understanding is

grad(central_fdm(5, 1), f, x)

requires x, an array specifying the scalar value spatial positions.

It is not clear to me as to what is the best way to set this up for a 3D scalar array.

Any advice on how to proceed would be appreciated.

2

It seems that you have a scalar array but not a scalar function defined?
If that is the case then check this post out.

3

Quite right. It is numerical data for which there is no simple scalar function.
Specifically, a 3D image volume consisting of a stack of 2D radiological [CT] image slices.

Thanks for your suggestion and link.
However, my implementation, as below, takes far too long to be of use.

#=
rows_coords: Vector{Float64} (512,)
cols_coords: Vector{Float64} (512,)
slices_coords: Vector{Float64} (802,)
=#
nodes = (rows_coords, cols_coords, slices_coords)

#=
seriesImgVlm: Array{Float64, 3}(512, 512, 802)
=#
interpolatedImgVlm = interpolate(nodes, seriesImgVlm, Gridded(Linear()))
#=
n_rows:  512
n_cols: 512
n_slices: 802
=#
gradientMgntdSqrdImgVlm = zeros(Float64, n_rows, n_cols, n_slices)

@time begin
    for k = 1:n_slices
        slice_coord = slices_coords[k]

        for j = 1:n_cols
            col_coord = cols_coords[j]

            for i = 1:n_rows
                row_coord = rows_coords[i]
                gradientImgVlm = 
                    Interpolations.gradient(interpolatedImgVlm, row_coord, col_coord, slice_coord)
                gradientMgntdSqrdImgVlm[i, j, k] = dot(gradientImgVlm, gradientImgVlm)
            end
        end
    end
end

The timer reports, for the 2nd run - after compilation is done:

134.354073 seconds (1.47 G allocations: 43.902 GiB, 6.24% gc time)

As I am still learning Julia, coming from a C++/Python background, I would appreciate any advice on how to improve the performance of calculating the magnitude of the gradient of a 3D array. The code I am rewriting from C++ requires repeated gradient calculations, so speed is essential. I had first considered using diff(), but would like to avoid having to pad the arrays to deal with the boundaries.
Also, the use of Interpolations.jl or FiniteDifferences.jl was recommended in Julia discourse.

4

First and foremost, it’s important to make sure you’re benchmarking within a function to avoid slowdowns associated with global variables.

If your grid points are uniformly- but not equally-spaced, you’ll want to use scaled b-splines, which take advantage of the uniformity along each axis.

Taking these points into account, here’s an implementation of what I think you’re trying to do:

function abs²grad(img, grid)
    itp = scale(interpolate(img, BSpline(Linear())), grid...)
    map(Iterators.product(grid...)) do (x, y, z)
        ∇img = Interpolations.gradient(itp, x, y, z)
        dot(∇img, ∇img)
    end
end

julia> img = rand(512, 512, 802);

julia> grid = LinRange.(0, rand(), size(img));

julia> @time abs²grad(img, grid);
 10.400567 seconds (4 allocations: 3.133 GiB, 0.92% gc time)
1 Like
5

Thanks for your Julia link and code examples.

6

Our new package EquivariantOperators.jl does exactly this:

using EquivariantOperators
cell = [dx 0; 0 dy]
▽ = Del(cell)
▽(a)

You can also take div, curl, laplacian, custom green’s functions…

Tutorials on colab
Paper Preprint
Docs