clear;clc;
% part1 of hw #3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load the data
load p1.mat;
I1 = p1;
load p2.mat;
I2 = p2;

% compute centroid of each image
L1 = ones(size(I1));    % define regions
L1(find(I1==0)) = 0;    % assing 0 to non-image regions
stat1 = regionprops(L1, I1,'centroid');
x1 = stat1.Centroid(1); y1 = stat1.Centroid(2);
L2 = ones(size(I2));    % define regions
L2(find(I2==0)) = 0;    % assing 0 to non-image regions
stat2 = regionprops(L2, I2,'centroid');
x2 = stat2.Centroid(1); y2 = stat2.Centroid(2);

% print image centroid coordinates
fprintf(1, 'CI1 = (%.2f, %.2f)\n', x1, y1);
fprintf(1, 'CI2 = (%.2f, %.2f)\n', x2, y2);

% draw the centroids on the images
colormap(gray);
hold on; imagesc(I1); colormap(gray); plot(x1, y1, 'r*', 'MarkerSize', 8);hold off;
figure; hold on; imagesc(I2); colormap(gray); plot(x2, y2, 'r*', 'MarkerSize', 8);hold off;

% estimate the transformation between two images, where
%   T(x) = x + T where x belongs to I1
myT = I1 - I2;
% now plot the transformation
figure; hold on; imagesc(myT); colormap(gray); hold off;

% part2 of hw #3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% compute covariances, where cov(p) = E_p((x-c) * (x-c)')
% fistly apply centroid, i.e. shift elements
I1_C = zeros(size(I1));   % I1_C stands for I1 - C
I2_C = zeros(size(I2));
for i = 162:size(I1,1)
    for j = 173:size(I1,2)
        I1_C(i,j) = I1(i-161, j-172);
    end
end
for i = 173:size(I2,1)
    for j = 163:size(I2,2)
        I2_C(i,j) = I2(i-172, j-162);
    end
end

% cov should be 2 by 2
cov1 = I1_C * I1_C';
cov2 = I2_C * I2_C';


% part3 of hw #3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% compute principal axis and rotation angle
[mySigma1 myD1] = eig(cov1);    % mySigma's are orthogonal rotation matrices
[mySigma2 myD2] = eig(cov2);    % D's are eigenvalues
myR = mySigma1 * mySigma2';
myD = acos(myR);


% part4 of hw #3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% use the translation and rotation parameters to obtain image p3
% by transforming image p2 into coordinate system of p1
myP3 = mySigma1 * mySigma2' * I2_C
for i = 162:size(I1,1)
    for j = 173:size(I1,2)
        I1_C(i,j) = I1(i-161, j-172);
    end
end