MATLAB: Manually determining the gradient of a surface

beginnerconvolutiongradientkernelmanual

Hi , I’m trying to write a function that can manually determine the gradient of a surface (represented by a multivariable function that I have created) at a particular point (x,y)
The code I've written is shown below.
function [pdx,pdy] = comp_gradient(x,y)
kernelx = [2.5,0,-2.5];
kernely = [2.5;0;-2.5];
z = cos(x/2).* cos(y) + (y/10) – (x/5);
pdx = conv2(kernelx,z,'same');
pdy = conv2(kernely,z,'same');
end
I'm not sure where to go from here, my output is clearly incorrect.
Any help would be greatly appreciated.

Best Answer

Your attempt as you posted it is not actually that far off.
First, we need to know what the spacing is between elements, and in each direction. So what is an estimator of the first derivative of a function? Assume we have a function defined on a uniform spacing, with spacing dx and dy in the x and y directions.
First create a grid of points. Note that I used a different spacing in x and y to make things clear.
x = linspace(-1,1,20);
dx = x(2) - x(1);
y = linspace(-2,2,20);
dy = y(2) - y(1);
[X,Y] = meshgrid(x,y);
% and a function
fun = @(x,y) cos(x/2).* cos(y) + (y/10) - (x/5);
Z = fun(X,Y);
Now, you need to remember that meshgrid variex x as the SECOND index. I used meshgrid in this because we tend to think of X as the horizontal coordinate, Y as the vertical. But we will need to be careful, as the horizontal coordinate of an array is the second one. Anyway...
What is the first derivative in the x direction, or an estimate thereof? You used central differences, so lets stick with that. We would have something like this:
df/dx = (F(y, x + dx) - F(y , x - dx))/(2*dx)
df/dy = (F(y + dy, x) - F(y - dy , x))/(2*dy)
We can achieve that using a convolution, as you tried. The convolution kernel for the x and y gradients would be
kernelx = [-1 0 1]/(2*dx);
kernely = [-1;0;1]/(2*dy);
Now, how to use conv? Be careful, as edge effects will cause problems. conv2 will insert zeros for those terms that fall off the edge of the world. So I would strongly recommend the use of 'valid'. 
pdx = conv2(Z,kernelx,'valid');
pdy = conv2(Z,kernely,'valid');
Note that this will cause the partial derivative estimate arrays to be 2 elements smaller in the indicated directions. We see that in the respective sizes of pdx and pdy.
size(pdx)
ans =
    20    18
size(pdy)
ans =
    18    20
If you wish, you could patch the edges, perhaps inserting a lower order approximation as new first and last columns for x, and new first and last rows for y.
pdx = [conv2(Z(:,1:2),[-1 1]/dx,'valid'), pdx, conv2(Z(:,[end-1,end]),[-1 1]/dx,'valid')];
size(pdx)
ans =
    20    20
surf(X,Y,pdx)
Then do something similar for the y partial.
So, as I said, you really were reasonably close in how you started.
Related Question