function [corners,edges,imcenter,I]=vharris(I_in,t1,t2,circsize)

% [corner_list,corners,lambda2,edges,both]=vharris(I,t1,t2)
% based on Harris & Stephes (1988)
% computa bordas e quinas a partir de autovalores lambda1, lambda2 e
% autovetor correspondente a lambda2.
% inputs: I imagem, t1 limiar
% para deteccao de quinas, t2 limiar de detccao de bordas
% outputs: corner_list: lista de coordenadas de quinas, corners:matriz de quinas, edges:matriz de arestas, both: matriz
% da combinacao das duas anteriores
% implementado por Lucas Pinto

% get inputs
if size(I_in,3)==3
    cI=double(I_in);
    I_in=double(rgb2gray(I_in));
else I_in=double(I_in);
end

if nargin==2
    t1=25;
    t2=5;
elseif nargin==1
    t1=25;
    t2=5;
end

% convolve dX^2, dY^2, dXY with gaussian kernel of varying size given by
% the relations:
% RF size(deg)=0.57*excentricity^0.69 (galago, Rosa et al. 1997)
% 1 pixel=0.25 degrees, (480 pixels=120 deg, size of human binocular field)
% minimum allow N=2 (window size=2*N+1=5), to avoid gaussian aliasing
% which yields: N=round((((0.57*excentricity^0.69)/4)/2)-1+3);

%apply circular mask
[I,imcenter]=circmask(I_in,circsize);

% calculate image gradient
[dX,dY]=gradient(I);

%calculate size relations
pixel2degree=circsize/120;
max_excentricity=round(sqrt((1-imcenter(1,1))^2+(1-imcenter(1,2))^2));
max_N=round((((0.57*max_excentricity^0.69)/pixel2degree)/2)-1+3);

%dX^2
x2=dX.^2;
dX2=zeros(size(x2));
Ir=zeros(size(I));

for i=1+max_N:size(I,1)-max_N
	for j=1+max_N:size(I,2)-max_N
        excentricity=round(sqrt((i-imcenter(1,1))^2+(j-imcenter(1,2))^2));
        N=round((((0.57*excentricity^0.69)/pixel2degree)/2)-1+3);
        sigma=(2*N+1)/5;
        Gr=exp(-(((-N:N).^2)./(2*sigma^2)));
        Gr=Gr./sum(Gr);
		for k=-N:N
			Ir(i,j)=Ir(i,j)+(Gr(k+N+1)*x2(i,j+k));
		end
	end
end

for i=1+max_N:size(I,1)-max_N
	for j=1+max_N:size(I,2)-max_N
        excentricity=round(sqrt((i-imcenter(1,1))^2+(j-imcenter(1,2))^2));
        N=round((((0.57*excentricity^0.69)/pixel2degree)/2)-1+3);
        sigma=(2*N+1)/5;
        Gc=exp(-(((-N:N).^2)./(2*sigma^2)));
        Gc=(Gc./sum(Gc))';
		for h=-N:N
			dX2(i,j)=dX2(i,j)+(Gc(h+N+1)*Ir(i+h,j));
		end
	end
end

%dY^2
y2=dY.^2;
dY2=zeros(size(y2));
Ir=zeros(size(I));

for i=1+max_N:size(I,1)-max_N
	for j=1+max_N:size(I,2)-max_N
        excentricity=round(sqrt((i-imcenter(1,1))^2+(j-imcenter(1,2))^2));
        N=round((((0.57*excentricity^0.69)/pixel2degree)/2)-1+3);
        sigma=(2*N+1)/5;
        Gr=exp(-(((-N:N).^2)./(2*sigma^2)));
        Gr=Gr./sum(Gr);
		for k=-N:N
			Ir(i,j)=Ir(i,j)+(Gr(k+N+1)*y2(i,j+k));
		end
	end
end

for i=1+max_N:size(I,1)-max_N
	for j=1+max_N:size(I,2)-max_N
        excentricity=round(sqrt((i-imcenter(1,1))^2+(j-imcenter(1,2))^2));
        N=round((((0.57*excentricity^0.69)/pixel2degree)/2)-1+3);
        sigma=(2*N+1)/5;
        Gc=exp(-(((-N:N).^2)./(2*sigma^2)));
        Gc=(Gc./sum(Gc))';
		for h=-N:N
			dY2(i,j)=dY2(i,j)+(Gc(h+N+1)*Ir(i+h,j));
		end
	end
end

%dXdY
xy=dX.*dY;
dXY=zeros(size(xy));
Ir=zeros(size(I));

for i=1+max_N:size(I,1)-max_N
	for j=1+max_N:size(I,2)-max_N
        excentricity=round(sqrt((i-imcenter(1,1))^2+(j-imcenter(1,2))^2));
        N=round((((0.57*excentricity^0.69)/pixel2degree)/2)-1+3);
        sigma=(2*N+1)/5;
        Gr=exp(-(((-N:N).^2)./(2*sigma^2)));
        Gr=Gr./sum(Gr);
		for k=-N:N
			Ir(i,j)=Ir(i,j)+(Gr(k+N+1)*xy(i,j+k));
		end
	end
end

for i=1+max_N:size(I,1)-max_N
	for j=1+max_N:size(I,2)-max_N
        excentricity=round(sqrt((i-imcenter(1,1))^2+(j-imcenter(1,2))^2));
        N=round((((0.57*excentricity^0.69)/pixel2degree)/2)-1+3);
        sigma=(2*N+1)/5;
        Gc=exp(-(((-N:N).^2)./(2*sigma^2)));
        Gc=(Gc./sum(Gc))';
		for h=-N:N
			dXY(i,j)=dXY(i,j)+(Gc(h+N+1)*Ir(i+h,j));
		end
	end
end

% eigenvalues and 1st eigenvector from matrix trace and determinant
trace=dX2+dY2;
det=(dX2.*dY2)-dXY.^2;
lambda1=0.5*(trace+ sqrt(trace.^2-4*det));
lambda2=0.5*(trace-sqrt(trace.^2-4*det));
v1=zeros(size(I,1),size(I,2),2);
v1(:,:,1)=dX2-lambda2;
v1(:,:,2)=dXY;
norm=sqrt(v1(:,:,1).^2+v1(:,:,2).^2);
norm(find(norm==0))=1;
v1(:,:,1)=v1(:,:,1)./norm;
v1(:,:,2)=v1(:,:,2)./norm;

% detect edges and corners
Es=zeros(size(I));
Eo=zeros(size(I));
corners=zeros(size(I));
edges=zeros(size(I));
lambda1=lambda1.*255/(max(max(lambda1)));
lambda2=lambda2.*255/(max(max(lambda2)));

for i=1:size(I,1)
    for j=1:size(I,2)
        if lambda2(i,j)<t1
            Es(i,j)=lambda1(i,j);
            Eo(i,j)=rad2deg(real(atan(v1(i,j,2))/(v1(i,j,1)+0.00001)));
            if Eo(i,j)>360|Eo(i,j)<-360
                Eo(i,j)=rem(Eo(i,j),360);
            end
            if Eo(i,j)>180 & Eo(i,j)<360
                Eo(i,j)=Eo(i,j)-360;
            end
        elseif lambda2(i,j)>=t1
                corners(i,j)=1;
        end
    end
end

% corner suppression
[k,w]=find(corners);
initial_corner_list=[k w];

for m=1:size(k,1)
    lambda2_list(m)=lambda2(k(m),w(m));
end
lambda2_list=lambda2_list';

[sorted_lambda2,sorted_i]=sort(lambda2_list,'descend');
sorted_corners=zeros(size(initial_corner_list));

for n=1:size(sorted_corners)
    sorted_corners(n,1:2)=initial_corner_list(sorted_i(n),1:2);
end

for p=1:size(sorted_corners,1)
    for q=1:size(sorted_corners,1)-p
        if (sorted_corners(p+q,1)>(sorted_corners(p,1)-N) & sorted_corners(p+q,1)<(sorted_corners(p,1)+N)) & (sorted_corners(p+q,2)>(sorted_corners(p,2)-N) & sorted_corners(p+q,2)<(sorted_corners(p,2)+N))  
        sorted_corners(p+q,1:2)=0;
        end
    end
end

[k,w]=find(sorted_corners);
corner_list=zeros(size(k,1)/2,2);
for i=1:size(k,1)/2
    corner_list(i,1:2)=sorted_corners(k(i),1:2);
end

corners=zeros(size(corners));
for i=1:size(corner_list)
    corners(corner_list(i,1),corner_list(i,2))=1;
end

% edge nonmaximum suppression & thresholding 

for i=2:size(Es,1)-1
    for j=2:size(Es,2)-1
        if Eo(i,j)<=22.5 & Eo(i,j)>=-22.5 | Eo(i,j)<=-157.5 | Eo(i,j)>=157.5
            if Es(i,j)<Es(i,j-1) | Es(i,j)<Es(i,j+1) | Es(i,j)<t2
                edges(i,j)=0;
            else edges(i,j)=1;
            end
        elseif Eo(i,j)>22.5 & Eo(i,j)<67.5 | Eo(i,j)>-157.5 & Eo(i,j)<-112.5
            if Es(i,j)<Es(i+1,j-1) | Es(i,j)<Es(i-1,j+1) | Es(i,j)<t2
                edges(i,j)=0;
            else edges(i,j)=1;
            end
        elseif Eo(i,j)>=67.5 & Eo(i,j)<=112.5 | Eo(i,j)>=-112.5 & Eo(i,j)<=-67.5
            if Es(i,j)<Es(i-1,j) | Es(i,j)<Es(i+1,j) | Es(i,j)<t2 
                edges(i,j)=0;
            else edges(i,j)=1;
            end
        elseif Eo(i,j)>112.5 & Eo(i,j)<=157.5 | Eo(i,j)>-67.5 & Eo(i,j)<=-22.5
            if Es(i,j)<Es(i-1,j-1) | Es(i,j)<Es(i+1,j+1) | Es(i,j)<t2
                edges(i,j)=0;
            else edges(i,j)=1;
            end
        end
    end
end

figure;
imshow(logical(corners));
figure;
imshow(logical(edges));

%plot together
%ccorners=zeros(size(corners,1),size(corners,2),3);
%cedges=zeros(size(edges,1),size(edges,2),3);
%[k,w]=find(corners);
%for i=1:size(k)
%ccorners(k(i),w(i),1:3)=[255 0 0];
%end
%[k,w]=find(edges);
%for i=1:size(k)
%cedges(k(i),w(i),1:3)=[0 0 255];
%end
%total=uint8(cI+ccorners+cedges);
%figure;
%imshow(total)