%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Taken from P.Kovesi Functions   HARRIS MEASURE
%

close all
clear

SIGMA_gaussian=4

y=im2double(imread('image_lena512.png'));
figure(1), imshow(y),title('original image');

dx = [-1 0 1; -1 0 1; -1 0 1];   % Derivative masks
dy = dx';

Ix = conv2(y, dx, 'same');      % Image derivatives
Iy = conv2(y, dy, 'same');

g = fspecial('gaussian',max(1,fix(3*SIGMA_gaussian)+1), SIGMA_gaussian);

Ix2 = conv2(Ix.^2, g, 'same'); % Smoothed squared image derivatives
Iy2 = conv2(Iy.^2, g, 'same');
Ixy = conv2(Ix.*Iy, g, 'same');

cim = (Ix2.*Iy2 - Ixy.^2)./(Ix2 + Iy2 + eps);

%%  remove boundaries
BORDER=10
cim(1:BORDER,:)=0;
cim(end-BORDER:end,:)=0;
cim(:,end-BORDER:end)=0;
cim(:,1:BORDER)=0;

figure(2), imshow(cim,[]),title('Harris measure');

%% Thresholding (possibly)

T=mean(cim(:));
CIM=cim;
CIM(find(cim<T))=0;

figure(3), imshow(CIM,[]),title('Harris measure');

%% perform nonlocal maximum suppression
support=true(11);
maxima=ordfilt2(CIM,sum(support(:)),support);
[loc_x,loc_y]=find((cim==maxima).*(CIM>0));

figure(3), imshow(y,[]), hold on, plot(loc_y,loc_x,'r+')