function rotation_invariance

img = imread("patch.png");

% circle mask
[h,w] = size(img);
s = min(h,w);
radius = s/2;

mask = zeros(h,w);

for y=1:h
    for x=1:w
        dx = x - w/2;
        dy = y - h/2;
        r2 = dx*dx + dy*dy;

        if r2 < radius*radius
            mask(y,x) = 1;
        end
    end
end

sobel = fspecial('sobel');

k=0;
for i=0:45:360-45
    img2 = imrotate(img, i, 'bilinear', 'crop');

    dx = conv2(img2, sobel', 'same');
    dy = conv2(img2, sobel, 'same');

    dx = dx .* mask;
    dy = dy .* mask;

    A = [dx(:) dy(:)];
    A = A' * A;

    % to prevent numerical overflow, the values can get VERY large
    A /= sum(mask(:));

    % calculate eigenvalue/eigenvector algebraically
    t = trace(A);
    d = det(A);
    eig1 = t*0.5 + sqrt(t*t*0.25 - d);

    angle = atan2(A(2,1), eig1 - A(2,2)) * 180 / pi;

    img_r1 = imrotate(img2, angle, 'bilinear', 'crop');
    img_r2 = imrotate(img2, angle + 180, 'bilinear', 'crop');

    k++;

    subplot(3,8,k);
    imshow(img2);

    subplot(3,8,k + 8);
    imshow(img_r1);
    title(sprintf("%.0f", angle));

    subplot(3,8,k + 16);
    imshow(img_r2);
    title(sprintf("%.0f", angle + 180));
end
