0

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. enter image description here

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.

MathIsHard
  • 2,733
  • 16
  • 47
  • 1
    What is the question? Are you looking for a suggestion on how to improve this? Are you looking for comments on what you have done? One thing is that you should have Cxx.*Cyy-(Cxy).^2. You can also use Matlab's "gradient" to compute the Hessian efficiently: [Cx, Cy] = gradient(Image); [Cxx, Cxy] = gradient(Cx); [Cyx, Cyy] = gradient(Cy); (Note that Cxy = Cyx). – Jeremy Upsal May 03 '18 at 20:57
  • I guess I am looking for some code that will work. I have been working on this for a while. I need this further some research I am trying to do... – MathIsHard May 03 '18 at 21:10
  • @JeremyUpsal Thank you. I did have that term squared. I seem to have lost it when I was copying the code over... Thank you for pointing that out for me. – MathIsHard May 03 '18 at 21:12
  • Another thing that needs to change is that you need Cx = Cy = 0 before using the second derivative test. So you should probably search for points first (using find) and then apply the second derivative test. Of course it's unlikely that it will be exactly 0, so you probably need some tolerance. I see that you were using Cx*Cy<0 before, why did you give that up? – Jeremy Upsal May 03 '18 at 21:17
  • I gave it up because I felt like I should be using the test itself... I think maybe I am misunderstanding how the test works. Are you saying that I want to test for Cx=Cy=0 (with some tolerance) and then apply the test to those points? – MathIsHard May 03 '18 at 21:21

1 Answers1

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. enter image description here

Nakini
  • 128
  • What I am seeing is that it picks up max and mins not close to the boundary. Consider this Image=0.1*[1 1 3 1 1;1 4 2 4 2; 1 2 1 2 1; 1 2 1 3 2; 1 2 1 3 2]; It is similar to the test image but I added another row moving the max to be more interior, and it picks up the max value now . – MathIsHard May 04 '18 at 16:26
  • @MathIsHard, the problem is that the code does not consider numerical approximations. Some critical points which are negative but close to zero should further be examined whether they are saddle points or not. – Nakini May 04 '18 at 18:41
  • Ah I see. Thank you. I will look into that. :) – MathIsHard May 04 '18 at 18:47
  • do I need to check first for the critical points and then run the test just on those points or do I test all the points? I’m a bit confused on that point. – MathIsHard May 04 '18 at 18:49
  • 2
    Isn't it possible to pick up values that are not true saddle points with this test? (i.e. points at which Cx and Cy aren't both 0?) – Jeremy Upsal May 04 '18 at 20:04
  • I would think so... But I am still getting max and min included with the test. I am trying to figure out how fix the problem. I tried to only test points that have Cx=Cy=0 approx. but I isn't working as of yet. – MathIsHard May 05 '18 at 03:18