%%
clear; load('visiblehuman.mat')
%%
mri = double(head_mri);
frozen = double(head_frozen);

[X, Y] = meshgrid(1:size(mri,2),1:size(mri,1));
mri_sum = sum(mri(:));
centerOfMassX_MR = sum(mri(:) .* X(:)) / mri_sum;
centerOfMassY_MR = sum(mri(:) .* Y(:)) / mri_sum;

mu_mri = [centerOfMassX_MR, centerOfMassY_MR]

[X, Y] = meshgrid(1:size(frozen,2),1:size(frozen,1));
frozen_sum = sum(frozen(:));
centerOfMassX_frozen = sum(frozen(:) .* X(:)) / frozen_sum;
centerOfMassY_frozen = sum(frozen(:) .* Y(:)) / frozen_sum;

mu_frozen = [centerOfMassX_frozen, centerOfMassY_frozen]

%%
mri_sigma = zeros(2);
for i=1:size(mri,1)
    for j=1:size(mri,2)
        mri_sigma = mri_sigma + mri(i, j)*([i; j] - mu_mri')*([i; j] - mu_mri')';
    end
end
mri_sigma = mri_sigma / mri_sum

frozen_sigma = zeros(2);
for i=1:size(frozen,1)
    for j=1:size(frozen,2)
        frozen_sigma = frozen_sigma + frozen(i, j)*([i; j] - frozen_mri')*([i; j] - frozen_mri')';
    end
end
frozen_sigma = frozen_sigma / mri_sum