function [] = edge(file, sigma, nsize, T)

%Read image, convert to grayscale and then to double matrix
img = imread(file);
img2 = img;
gray = rgb2gray(img);
d = im2double(gray);

[row col] = size(d);

%Create the filter size based on sigma, make sure its size is odd
fsize = round(sigma*6);
if mod(fsize, 2) == 0 || fsize == 0
    fsize = fsize + 1;
  
half = floor(fsize / 2);    

%Centralize the Gaussian filter so point 0,0 is at center of filter
for x = -1*half:half
    for y = -1*half:half
        Gy(x+half+1,y+half+1) = dGaussian(x,sigma)*Gaussian(y,sigma);
        Gx(x+half+1,y+half+1) = dGaussian(y,sigma)*Gaussian(x,sigma);
    end
end

%Perform convolution
Fx = conv2(d,Gx,'same');
Fy = conv2(d,Gy,'same');

L = zeros(row,col,2);

for x = 1:row
    for y = 1:col
        F(x,y) = sqrt(Fx(x,y)^2 + Fy(x,y)^2);
        [x_l,x_h,y_l,y_h] = getNDim(x,y,row,col,nsize);
        Fx2 = 0;
        Fy2 = 0;
        Fxy = 0;
        
        %Compute the covariance matrix over neighborhood
        for i = x_l:x_h
            for j = y_l:y_h
                Fx2 = Fx2 + (Fx(i,j)^2);
                Fy2 = Fy2 + (Fy(i,j)^2);
                Fxy = Fxy + (Fx(i,j)*Fy(i,j));
            end
        end
        
        %Find smallest eigenvalue of covariance matrix
        C = [Fx2 Fxy; Fxy Fy2];
        E = eig(C);
        minEig = min(E);
        
        %Find all corner pixels
        if minEig > T
            L(x,y,1) = minEig;
            L(x,y,2) = 1;
            img2(x,y,1) = 255;
            img2(x,y,2) = 0;
            img2(x,y,3) = 0;
        end
    end
end

%Suppress neighborhood corner pixels whose min eigenvalue is less than the
%current pixel's min eigenvalue
for x = 1:row
    for y = 1:col
        [x_l,x_h,y_l,y_h] = getNDim(x,y,row,col,nsize);
        for i = x_l:x_h
            for j = y_l:y_h    
                if L(i,j,1) < L(x,y,1)
                    L(i,j,2) = 0;
                end
            end
        end
    end
end

%Draw box around center of nonsurpressed corner pixels
for x = 1:row
    for y = 1:col
        if L(x,y,2) == 1
            [x_l,x_h,y_l,y_h] = getNDim(x,y,row,col,nsize);
            for i = x_l:x_h
                j = y_l;
                img(i,j,1) = 255;
                img(i,j,2) = 255;
                img(i,j,3) = 0;
                j = y_h;
                img(i,j,1) = 255;
                img(i,j,2) = 255;
                img(i,j,3) = 0;                
            end
            for j = y_l+1:y_h-1
                i = x_l;
                img(i,j,1) = 255;
                img(i,j,2) = 255;
                img(i,j,3) = 0;
                i = x_h;
                img(i,j,1) = 255;
                img(i,j,2) = 255;
                img(i,j,3) = 0;                
            end
        end
    end
end

%Convert matrices to images
GradX = mat2gray(abs(Fx));
GradY = mat2gray(abs(Fy));
GradM = mat2gray(F);

%Save image steps
params = ['-corner-' num2str(sigma) '-' num2str(nsize) '-' num2str(T) '-'];
nfile = strrep(file, '.jpg', '');
imwrite(GradX,[nfile params 'fgx.jpg']);
imwrite(GradY,[nfile params 'fgy.jpg']);
imwrite(GradM,[nfile params 'fgm.jpg']);
imwrite(img2,[nfile params 'fc.jpg']);
imwrite(img,[nfile params 'nms.jpg']);

imshow(img)
end