How do I compute the gradient vector of pixels in an image? I'm trying to find the curvature of the features in an image and I was advised to calculate the gradient vector of pixels. So if the matrix below are the values from a grayscale image, how would I go about calculating the gradient vector for the pixel with the value '99'?
21 20 22 24 18 11 23
21 20 22 24 18 11 23
21 20 22 24 18 11 23
21 20 22 99 18 11 23
21 20 22 24 18 11 23
21 20 22 24 18 11 23
21 20 22 24 18 11 23

Apologies for asking such an open ended question, I've never done much maths and am not sure how to start tackling this.
 A: In Python you can use the numpy.gradient function to do this. This said function uses central differences for the computation, like so:
\begin{eqnarray}
\nabla_x I (i, j) = \frac{I(i + 1, j) - I(i - 1, j)}{2}, \hspace{.5em}\nabla_y I (i, j) = \frac{I(i, j+1) - I(i, j-1)}{2}.
\end{eqnarray}
Here is a code snippet for your specific image:
import numpy as np
import matplotlib.pyplot as plt

# load image
img = np.array([[21.0, 20.0, 22.0, 24.0, 18.0, 11.0, 23.0],
                [21.0, 20.0, 22.0, 24.0, 18.0, 11.0, 23.0],
                [21.0, 20.0, 22.0, 24.0, 18.0, 11.0, 23.0],
                [21.0, 20.0, 22.0, 99.0, 18.0, 11.0, 23.0],
                [21.0, 20.0, 22.0, 24.0, 18.0, 11.0, 23.0],
                [21.0, 20.0, 22.0, 24.0, 18.0, 11.0, 23.0],
                [21.0, 20.0, 22.0, 24.0, 18.0, 11.0, 23.0]])
print "image =", img

# compute gradient of image
gx, gy = np.gradient(img)
print "gx =", gx
print "gy =", gy

# plotting
plt.close("all")
plt.figure()
plt.suptitle("Image, and it gradient along each axis")
ax = plt.subplot("131")
ax.axis("off")
ax.imshow(img)
ax.set_title("image")

ax = plt.subplot("132")
ax.axis("off")
ax.imshow(gx)
ax.set_title("gx")

ax = plt.subplot("133")
ax.axis("off")
ax.imshow(gy)
ax.set_title("gy")
plt.show()


To answer your specific question, the gradient (via central differences!) of the image at pixel with value $99$ is $0$ along the $x$ axis and $-2$ along the $y$ axis.
A: Suppose the image is continuous and differentiable in $x$ and $y$. Then $I(x,y)$ is the value of the pixel at each $(x,y)$, i.e. $I: \mathbb{R}^2 \mapsto \mathbb{R}$. Recall that the gradient at a point $(u,v)$ is:
$$
\nabla I(u,v) = \begin{bmatrix} \frac{\partial I}{\partial x}(u,v)  \\ \frac{\partial I}{\partial y}(u,v) \end{bmatrix}
$$
Given a discrete grid, you should approximate the partial derivative in $x$ and $y$ directions using finite difference approximations at the point of interest.
Assume your function $I$ is sampled over points $\{1, \ldots, 7 \} \times \{1, \ldots, 7 \}$ in image-coordinates, i.e. $I(1,1) = 21$, $I(1,7) = 23$, etc... So you're looking for the gradient at $(4,4)$. If you assume the resolution between points is 1, then the forward difference approximation in the $x$ direction gives:
$$
\frac{\partial I}{\partial x}(4,4) \approx I(5,4) - I(4,4) = 24 - 99
$$
Do the same in $y$ to obtain the full gradient at the point.
A: The theoretical thing you may want to read up on is convolution and especially discrete convolution.

Since there already are practical approaches in a few computer languages I might as well add one that will work in Matlab & GNU Octave. We start out by feeding in our image. It is typed in as a matrix:
lImage = [21,20,22,24,18,11,23;
          21,20,22,24,18,11,23;
          21,20,22,24,18,11,23;
          21,20,22,99,18,11,23;
          21,20,22,24,18,11,23;
          21,20,22,24,18,11,23;
          21,20,22,24,18,11,23];

Then the command "conv2" creates a 2 dimensional convolution, first we can do one in the x direction by feeding it a row vector with the [1,0,-1] filter:
dx = conv2(lImage,[1,0,-1],'valid')

which gives us the output:
dx =

   1    4   -4  -13    5
   1    4   -4  -13    5
   1    4   -4  -13    5
   1   79   -4  -88    5
   1    4   -4  -13    5
   1    4   -4  -13    5
   1    4   -4  -13    5 

and the y direction we can easily do by just applying a transpose (denoted by ') on the x filter row vector:
dy = conv2(lImage,[1,0,-1]','valid')

gives us the output:
dy= 

   0    0    0    0    0    0    0
   0    0    0   75    0    0    0
   0    0    0    0    0    0    0
   0    0    0  -75    0    0    0
   0    0    0    0    0    0    0

