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

Specific Domains Numerics
1

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
]
2

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

3

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)
	grads = zeros(size(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
	return grads
end
4

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
5

I guess central difference is a more accurate name. Essentially numpy.gradient
https://numpy.org/doc/stable/reference/generated/numpy.gradient.html

6

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
7

They form a curve in 3D space

8

Something like this in Makie

image
image1270×880 97.3 KB

9

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
10

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[[1], :]; dv]
    a .+= [dv; dv[[end], :]]
    a
end
11

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