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 
myT = [x1 y1] - [x2 y2];

%% part2 of hw #3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all; % close figures of previous part
% compute covariances, where cov(p) = E_p((x-c) * (x-c)')
cov1 = zeros(2,2);  
cov2 = zeros(2,2);
for i = 1:size(I1,1)
    for j = 1:size(I1,2)
        cov1 = cov1 + I1(i,j) * ([(i - x1);(j - y1)] * [(i - x1) (j - y1)]);
    end
end
for i = 1:size(I2,1)
    for j = 1:size(I2,2)
        cov2 = cov2 + I2(i,j) * ([(i - x2);(j - y2)] * [(i - x2) (j - y2)]);
    end
end
cov1 = cov1 / (size(I1,1) * size(I1,2));
cov2 = cov2 / (size(I2,1) * size(I2,2));

%% 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

turnangle = 0.0774 *pi;
myTranslation = [1 0 0; 0 1 0; -x2 -y2 1];
myTranslation2 = [1 0 0; 0 1 0; x1 y1 1];
myRotation = [cos(turnangle) -sin(turnangle) 0; sin(turnangle) cos(turnangle) 0; 0 0 1];

singleTrans = myTranslation * myRotation' * myTranslation2;

T = maketform('affine', singleTrans);
I3 = imtransform(I2,T,'bilinear', 'XData', [1 size(I1,2)],'YData',[1 size(I1,1)]);
diffImage = abs(I3-I1);
figure;imagesc(diffImage); colormap(gray);




