Program to find Saddle points matlab

Question

I am trying to write a program that will take in a matrix (lattice) that defines a surface and find all the saddle points. I was trying to use finite differences and the second partial derivative test to find the saddle points. I have tried all the methods I have been able to find online without much luck.

This is my code so far.

Image = getSurface([0.5 0.5],1)
    %Image=[1 1 3; 4 2 4; 2 1 3]
    Image=[1 1 3 1;4 2 4 2;2 1 2 1; 2 1 3 2];
    h=0.05;
    Cxx=diff(Image,2,1)/h^2; 
    Cxx(end+2,:)=nan;
    Cyy=diff(Image,2,2)/h^2; 
    Cyy(:,end+2)=nan;
    Cx=diff(Image,1,1)/h;
    Cx(end+1,:)=nan;
    Cxy=diff(Cx,1,2)/h;
    Cxy(:,end+1)=nan;

    [i,j]=find(Cxx.*Cyy - Cxy^2 < 0);
    %[i,j]=find(Cxx.*Cyy< 0);

% hold on
% for r=1:1:length(i)
%     plot(Image(i(r)),Image(j(r)),'.')
%      plot(i(r),j(r),'*')
% end
% hold off

figure;
hold on;
surf(Image);
shading interp
for r=1:1:size(i,1)
    plot3(i(r),j(r),4, 'k*');
end
hold off;

Here is the output that I am getting so far. I seem to be getting one of the saddle points along with some of the max and min points.

I am trying to see if anyone has some code that will work for this? I am trying to further some research and this is the hurdle I am stuck at right now. Thank you.

Answer 1

First, as @Jeremy Upsal said, it's better to use gradient() to calculate the gradient and the hessian. Second, to discard all the point on the boundary, you have to make little changes in your find() function. Lastly, you don't need a for-loop to superimpose the saddle points on the image. Take a look at the following code:

Image=[1 1 3 1;4 2 4 2;2 1 2 1; 2 1 3 2];

% Evaluate gradient and hessian
[Cx, Cy] = gradient(Image);
[Cxx, Cxy] = gradient(Cx);
[Cyx, Cyy] = gradient(Cy);
D = Cxx.*Cyy - Cxy^2;

% Discard all the points on the boundary
[i,j]=find(D(2:end-1, 2:end-1) < 0);
i = i + 1; j = j + 1;

% Plot the surface and the saddle points
figure;
hold on;
surf(Image);
shading interp
plot3(i,j,repmat(4,[length(i),1]), 'k*');
hold off;

You will get the following figure.