image = load('boat.mat');
image = image.imageOrig;

image = mat2gray(image);

I = image;
figure, imshow(mat2gray(I)); %first image

dim = size(image);

sigma1 = 2;
sigma2 = 2;
k = 0.06;

gauss1 = fspecial('gaussian', [3 3], sigma1);
gauss2 = fspecial('gaussian', [3 3], sigma2);

image = imfilter(image, gauss1);

[Ix, Iy] = imgradientxy(image);

figure, imshow(Ix); %image gradient x-dirn
figure, imshow(Iy); %image gradient y-dirn

Ix2 = imfilter(Ix.^2, gauss2);
Iy2 = imfilter(Iy.^2, gauss2);
IxIy = imfilter(Ix.*Iy, gauss2);

E1 = zeros(dim);
E2 = zeros(dim);
R = zeros(dim);

for i = 2:dim(1)-1
    for j = 2:dim(2)-1
        M = [Ix2(i,j), IxIy(i,j); IxIy(i,j), Iy2(i,j)];
        e = eig(M);
        E1(i,j)=min(e);
        E2(i,j)=max(e);
        R(i,j) = ((Ix2(i,j)*Iy2(i,j)) - (IxIy(i,j))^2 ) - k*(Ix2(i,j) + Iy2(i,j)).^2;
    end
end

figure, imshow(mat2gray(E1));
figure, imshow(mat2gray(E2));

R = R*1000/(max(max(R)));

cornerMeasure = (R>10);

maxR = max(max(R));

figure, imshow(mat2gray(cornerMeasure) + I);