Two solutions that don't work because of MATLAB grader

clustering_method.m

target = imread('example_target.png');

% harris returns two objects: pts, cim. cluster them
[pts, cim] = harris(target, 1, 1500, true);
T = clusterdata(pts', 0.8);

% Find the 48 largest clusters
[n, bin] = hist(T, unique(T));
[~,idx] = sort(-n);
largest_clusters = bin(idx(1:48));

% plot cluster on image
figure
imshow(target);
hold on;
for i=largest_clusters'
    clus = pts(:,T==i);
    k = convhull(clus');
    plot(clus(1,k), clus(2,k), 'LineWidth', 2);
    hold on;
end

% Find the largest element in each cluster
figure
imshow(target);
hold on;
result = zeros(2,48);
for i=1:48
    cluster_index = largest_clusters(i);
    clus = pts(:,T==cluster_index);
    [~, max_index] = max(cim(sub2ind(size(cim), clus(2,:), clus(1,:))));
    result(:,i) = clus(:, max_index);
end 
scatter(result(1,:), result(2,:), 'r+');

% Sort result points. This will make the points correspond to the right
% order as long as the image isn't /too/ rotated
sort(result, 2)

imfill_method.m

I = imread('example_target.png');

threshold = 1500;

[pts, cim] = harris(I, 1, threshold, false);
binarized_cim = cim >= threshold;

% find connected regions by filling them
% need to store size of region and highest point
regions = zeros(0, 2);
for start_point=pts
    if binarized_cim(start_point(2), start_point(1))
        connected_pixels = bwselect(binarized_cim, start_point(1), start_point(2), 4);
        
        % calculate the size of the region
        area = sum(sum(connected_pixels));        
        
        % find the location of the max (in 1D)
        [~, argmax] = max(cim(connected_pixels));
        pixel_indices = find(connected_pixels, 1000);
        location = pixel_indices(argmax);
        
        % store the region. not efficient, but neither is matlab.
        regions = [regions; area location];
        
        % mask the region off
        binarized_cim = binarized_cim & ~connected_pixels;
    end
end

% get the largest regions, and save the best point in each to
% `maxima`.
sorted_regions = sortrows(regions);
maxima = sorted_regions(end-48+1:end,2);

% extract a point that's representative of each region
% sort the points to get an ordering. will work if rotated <45 deg
[x,y] = ind2sub(size(I), maxima);
sorted_points = [y,x];

% get the saddle points
WINDOW = 7
result = zeros(2,48);
for i=1:48
    x = sorted_points(i,1);
    y = sorted_points(i,2);
    result(:,i) = [x-WINDOW-1;y-WINDOW-1] + saddle_point(double(I(y-WINDOW:y+WINDOW,x-WINDOW:x+WINDOW)));
end

figure
imshow(target);
hold on;
% scatter(sorted_points(:,1), sorted_points(:,2), 'r+');
scatter(result(1,:), result(2,:), 'r+');

% ADD YOUR SADDLE_POINT CODE BELOW
function [pt] = saddle_point(I)
[height, width] = size(I);
xcoords = repmat(1:width, 1, height); % x coordinates
ycoords = repelem(1:height, 1, width); % y coordinates

X = [ xcoords.^2;
    xcoords.*ycoords;
    ycoords.^2;
    xcoords;
    ycoords;
    ones(1, height*width)]';
Y = I(sub2ind([height, width], ycoords, xcoords))';

% Least squares to find coefficients alpha through eta
coeff = pinv(X) * Y;

% Find intersection using equation in the paper
pt = -inv([2*coeff(1)  coeff(2);
            coeff(2)     2*coeff(3)]) * [coeff(4); coeff(5)];

%------------------
end