# What is the best way to calculate the numerical gradient of a matrix, representing multiple points in 3D space?

Given a matrix, where each row represents a point in 3D space and each column represents the x, y, and z location. What is the recommended way to calculate the gradient at each point (at each row)?

``````points = [
202  279   8
202  280   9
209  278  10
210  283  12
212  286  14
219  295  16
225  300  18
234  308  20
238  305  22
236  300  24
233  295  26
231  289  28
227  280  30
223  275  31
]
``````

I haven’t verified it, but I am assuming `diff()` works. If so, is this the recommended approach? Maybe `imgradients()` from Images.jl?

Here is a simple working example I came up with

``````function centraldiff(p1, p2, p3)
u = p3 - p2
v = p2 - p1
return (u - v) ./ (norm(p3 - p1))
end
``````
``````function gradients(points)
for i in axes(points, 1)
if i == 1
grads[i, :] = centraldiff(points[i, :], points[i, :], points[i+1, :])
elseif i == size(points, 1)
grads[i, :] = centraldiff(points[i-1, :], points[i, :], points[i, :])
else
grads[i, :] = centraldiff(points[i-1, :], points[i, :], points[i+1, :])
end
end
end
``````

I’m not sure if the concept of ‘gradient’ makes any sense in this context. If `points` is just a list of positions, what are you taking the grafient of?

The gradient should represent the rate of change of a function or field with respect to position. But you just have the positions, not the ‘thing’ of which you should find the gradient.

The rate of change of positions with respect to position should just be constant.

1 Like

I guess central difference is a more accurate name. Essentially `numpy.gradient`

That’s a scheme for estimating the gradient. But what are you taking the gradient of? I mean, where are the data whose sample points you have provided?

1 Like

They form a curve in 3D space

Gradients are for functions of multiple variables, but you have a 3D curve parameterized by a single variable.

You can use `diff` for this, sure.

2 Likes

If that’s the case then this is probably the simplest answer (modified from this)

``````function central_difference(v::AbstractMatrix)
dv = diff(v, dims=1) / 2
a = [dv[, :]; dv]
a .+= [dv; dv[[end], :]]
a
end
``````

That looks very inefficient, allocating multiple arrays where only one should be needed. (Sorry, it’s too late for me to suggest a fix now, from my phone, but I suggest an array comprehension, or pre-allocate + loop. This is one example where Matlab-style vectorization is very suboptimal.)

1 Like