clf;
clear all;

Iorig = rgb2gray(imread('Shakey.jpg'));
imshow(Iorig);
I = double(Iorig);

mfs = 16;

% Convolution Sobel
SobelX = [-1 0 1; -2 0 2; -1 0 1];
SobelY = [1 2 1; 0 0 0; -1 -2 -1];
Ix = conv2(SobelX,I);
Iy = conv2(SobelY,I);
subplot(2,2,1); imagesc(Ix); colormap('gray'); title('Ix','FontSize',mfs); axis equal
subplot(2,2,2); imagesc(Iy); colormap('gray'); title('Iy','FontSize',mfs); axis equal



% Computing Ix2, Iy2 and Ixy matrices
Ix2 = Ix.^2;
Iy2 = Iy.^2;
Ixy = Ix.*Iy;

G = fspecial('gaussian', [5 5], 1);

Sx2 = conv2(Ix2,G);
Sy2 = conv2(Iy2,G);
Sxy = conv2(Ixy,G);

k = 0.1; % Parametre kappa
R = Sx2.*Sy2-Sxy.^2-k*((Sx2+Sy2).^2);

% Remove the edges
Border = 10;
R = R(Border+1:end-Border,Border+1:end-Border);
subplot(4,3,7); imagesc(R>1e4); title('R>1e4','FontSize',mfs); 
subplot(4,3,8); imagesc(R>1e6); title('R>1e6','FontSize',mfs); 
subplot(4,3,9); imagesc(R>1e7); title('R>1e7','FontSize',mfs); 

% Trouver les coordonnees des points saillants
[CornerX CornerY] = find(R>10000000);

figure(4)
clf
imshow(Iorig);
hold on;
plot(CornerY+Border-3,CornerX+Border-3,'ro');









