# Q: FiniteDifferences | gradient of a 3D scalar array

**URL:** <https://discourse.julialang.org/t/q-finitedifferences-gradient-of-a-3d-scalar-array/64065>\
**Category:** Numerics\
**Tags:** question\
**Created:** [July 5, 2021, 7:22am UTC](https://discourse.julialang.org/t/q-finitedifferences-gradient-of-a-3d-scalar-array/64065 "2021-07-05T07:22:25Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Audrius-St](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/audrius-st/32/24175_2.png) [@Audrius-St](https://discourse.julialang.org/u/Audrius-St)\
**Post date:** [July 5, 2021, 7:22am UTC](https://discourse.julialang.org/t/q-finitedifferences-gradient-of-a-3d-scalar-array/64065/1 "2021-07-05T07:22:25Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [July 5, 2021, 11:06am UTC](https://discourse.julialang.org/t/q-finitedifferences-gradient-of-a-3d-scalar-array/64065/2 "2021-07-05T11:06:00Z")

</div>

It seems that you have a scalar array but not a scalar function defined?  
If that is the case then check [this post](https://discourse.julialang.org/t/differentiation-without-explicit-function-np-gradient/57784) out.

---

<div class="post-metadata">

**Author:** ![Audrius-St](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/audrius-st/32/24175_2.png) [@Audrius-St](https://discourse.julialang.org/u/Audrius-St)\
**Post date:** [July 12, 2021, 8:33pm UTC](https://discourse.julialang.org/t/q-finitedifferences-gradient-of-a-3d-scalar-array/64065/3 "2021-07-12T20:33:25Z")

</div>

> [@rafael.guerra](#):
>
> It seems that you have a scalar array but not a scalar function defined?

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.

> [@rafael.guerra](#):
>
> If that is the case then check [this post](https://discourse.julialang.org/t/differentiation-without-explicit-function-np-gradient/57784) out.

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

```julia
#=
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.

---

<div class="post-metadata">

**Author:** ![stillyslalom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stillyslalom/32/45687_2.png) [@stillyslalom](https://discourse.julialang.org/u/stillyslalom)\
**Post date:** [July 13, 2021, 12:58am UTC](https://discourse.julialang.org/t/q-finitedifferences-gradient-of-a-3d-scalar-array/64065/4 "2021-07-13T00:58:17Z")

</div>

First and foremost, it’s important to make sure you’re benchmarking within a function to avoid [slowdowns associated with global variables](https://docs.julialang.org/en/v1/manual/performance-tips/#Avoid-global-variables).

If your grid points are uniformly- but not equally-spaced, you’ll want to use [scaled b-splines](http://juliamath.github.io/Interpolations.jl/latest/control/#Scaled-BSplines-1), 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:

```julia
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-repl
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)

```

---

<div class="post-metadata">

**Author:** ![Audrius-St](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/audrius-st/32/24175_2.png) [@Audrius-St](https://discourse.julialang.org/u/Audrius-St)\
**Post date:** [July 14, 2021, 4:30pm UTC](https://discourse.julialang.org/t/q-finitedifferences-gradient-of-a-3d-scalar-array/64065/5 "2021-07-14T16:30:13Z")

</div>

Thanks for your Julia link and code examples.

---

<div class="post-metadata">

**Author:** ![pxshen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pxshen/32/38776_2.png) [@pxshen](https://discourse.julialang.org/u/pxshen)\
**Post date:** [August 13, 2022, 2:05am UTC](https://discourse.julialang.org/t/q-finitedifferences-gradient-of-a-3d-scalar-array/64065/6 "2022-08-13T02:05:53Z")

</div>

Our new package `EquivariantOperators.jl` does exactly this:

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

```

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

[Tutorials on colab](https://colab.research.google.com/drive/17JZEdK6aALxvn0JPBJEHGeK2nO1hPnhQ#scrollTo=4utBHhoeHzRY)  
[Paper Preprint](http://arxiv.org/abs/2108.09541)  
[Docs](https://aced-differentiate.github.io/EquivariantOperators.jl/)