function X = find_corners(I, sigma, n, t)

[F, D, Fx, Fy] = filtered_gradient(I, sigma);

%Compute the values to be used for covariance
Fx2 = Fx.^2;
Fxy = Fx.*Fy;
Fy2 = Fy.^2;

%Use convolutions to sum up the neighborhoods
SFx2 = conv2(Fx2, ones(n*2+1), 'same');
SFxy = conv2(Fxy, ones(n*2+1), 'same');
SFy2 = conv2(Fy2, ones(n*2+1), 'same');

%Compute the minimum eigenvalue for the cov matrix of each pixel
L = zeros(size(I));
for i=n:size(I,1)-n;
    for j=n:size(I,2)-n;
        C = [SFx2(i,j) SFxy(i,j);
             SFxy(i,j) SFy2(i,j) ];
        L(i,j) = min(eig(C));
    end;
end;

%Threshold the results
L(L < t) = 0;

%Display the min0eigenvalue matrix
figure;
colormap(gray);
imagesc(L);
drawnow;

tic;

%Limit the searched region to those with a valid neighborhood
n=n+2;
[idxI idxJ] = find(L > 0);
idxI2 = idxI( idxI >= n+1 & idxI <= size(I,1)-n-1 & idxJ >= n+1 & idxJ <= size(I,2)-n-1 );
idxJ2 = idxJ( idxI >= n+1 & idxI <= size(I,1)-n-1 & idxJ >= n+1 & idxJ <= size(I,2)-n-1 );

%Perform the non-maximum test
M = zeros(size(I));
for k=1:length(idxI2);
    %Get the index of the next item, and the neighboring range
    i = idxI2(k);
    j = idxJ2(k);
    
    rangeI = i-n:i+n;
    rangeJ = j-n:j+n;   
    
    %determine if this has the maximum value in the neighboring pixels
    if( L(i,j) > 0 && all(all(repmat(L(i,j), 2*n+1, 2*n+1) >= L(rangeI,rangeJ))) );
        M(i,j) = 1;
    end;
end;

toc;

%Display the corner points
figure;
colormap(gray);
imagesc(M);
drawnow;

X = M;