function feature_list = corner_detector (A_img,n) 
% function feature_list = corner_detector (image,n) 
% image as obtained by im=imread('filename');
% n is the desired number of features

if nargin<2 
    n=1000;
end

percentage_overall=20/100;
s=size(A_img);
pix_per_region=s(1)*s(2)/36;

n_hor_region=round(s(2)/sqrt(pix_per_region));
n_vert_region=round(s(1)/sqrt(pix_per_region));

region_w=s(2)/n_hor_region;
region_h=s(1)/n_vert_region;

features_per_region=round(n*(1-percentage_overall)/36);
strongest_features_overall = n - 36*features_per_region;

min_threshold_in_percentage = 1/10; % to be a feature, its value has to be over this percentage among all features in the image
show_resulting_img = 0; % 1: show the resulting img; 0 : don't show

%input image
%A_img = imread(img_file);
img_h = size(A_img,1);
img_w = size(A_img,2);
A  = double(rgb2gray(A_img));

%parameters for the gaussian
gaussian_size = 7;
sigma = .7;

%feature window calculation
del_A_1 = conv2(A,[-1,1],'same') ;
del_A_2 = conv2(A,[-1;1],'same') ;
del_A_1_1 = del_A_1 .* del_A_1;
del_A_2_2 = del_A_2 .* del_A_2;
del_A_1_2 = del_A_1 .* del_A_2;
gaussian = fspecial('gaussian',gaussian_size,sigma);
matrix_1_1 = conv2(del_A_1_1,gaussian,'same');
matrix_2_2 = conv2(del_A_2_2,gaussian,'same');
matrix_1_2 = conv2(del_A_1_2,gaussian,'same');
to_i = img_h+1-gaussian_size;
to_j = img_w+1-gaussian_size;
k = 0.04;
R = zeros(img_h,img_w);
length_L = to_i*to_j;
L = zeros(length_L,3);
ind_L = 1;
for i= 1:to_i
    for j=1:to_j
        M = [matrix_1_1(i,j),matrix_1_2(i,j);
            matrix_1_2(i,j),matrix_2_2(i,j)];
        tM = trace(M);
        dM = det(M);
        R(i,j) = dM-k*tM^2; %min(abs(eig(M)));
%        L(ind_L,1) = R(i,j); %min(abs(eig(M)));
%        L(ind_L,2) = round(i);
%        L(ind_L,3) = round(j);
%        ind_L = ind_L + 1;
    end
end
for i= 2:to_i-1
    for j=2:to_j-1
%        R(i,j) = dM-k*tM^2; %min(abs(eig(M)));
        if (R(i,j)==max(max(R(i-1:i+1,j-1:j+1))))
            L(ind_L,1) = R(i,j); %min(abs(eig(M)));
            L(ind_L,2) = round(i);
            L(ind_L,3) = round(j);
            ind_L = ind_L + 1;
        end
    end
end
ind_L

%find features
feature_index=1;
is_not_a_feature = ones(img_h,img_w);
[sorted_L_1, sorted_ind] = sort(L,1);
min_threshold = L(sorted_ind(ceil(length_L*min_threshold_in_percentage)),1);
region = zeros(ceil(img_h/region_h), ceil(img_w/region_w));
total_num_region = ceil(img_h/region_h)*ceil(img_w/region_w)*features_per_region;
counter_total_region = 0;
for i = 1:length_L
    elem = sorted_ind(length_L-i+1);
    m = L(elem,2);
    n = L(elem,3);
    if i <= strongest_features_overall
        is_not_a_feature(m,n) = 0;
        feature_list(feature_index,1:2)= [m n];
        feature_index = feature_index + 1;
    else 
        p = ceil(m/region_h);
        q = ceil(n/region_w);
        if region(p,q) < features_per_region
            region(p,q) = region(p,q) +1;
            if L(elem,1) > min_threshold
                is_not_a_feature(m,n) = 0;
                feature_list(feature_index,1:2)= [m n];
                feature_index = feature_index + 1;
            end    
            counter_total_region = counter_total_region +1;
            if counter_total_region >= total_num_region
                break;
            end
        end
    end
end
if show_resulting_img == 1
    hA = imshow(uint8(A.*is_not_a_feature));
end
return;