k(1:4) = 0;
[x, y] = meshgrid(1:size(Ix,2), 1:size(Ix,1));

%%% Initialize the values in the A matrix
A = zeros(4);
A(1,1) = sum(sum((x.*Ix+y.*Iy).^2));
A(1,2) = sum(sum((x.*Ix+y.*Iy).*(y.*Ix-x.*Iy)));
A(2,1) = A(1,2);
A(1,3) = sum(sum(Ix.*(x.*Ix+y.*Iy)));
A(3,1) = A(1,3);
A(1,4) = sum(sum(Iy.*(x.*Ix+y.*Iy)));
A(4,1) = A(1,4);
A(2,2) = sum(sum((y.*Ix-x.*Iy).^2));
A(2,3) = sum(sum(Ix.*(y.*Ix-x.*Iy)));
A(3,2) = A(2,3);
A(2,4) = sum(sum(Iy.*(y.*Ix-x.*Iy)));
A(4,2) = A(2,4);
A(3,3) = sum(sum(Ix.^2));
A(3,4) = sum(sum(Ix.*Iy));
A(4,3) = A(3,4);
A(4,4) = sum(sum(Iy.^2));

%%% Initialize the values in the B matrix
B = zeros(4,1);
B(1,1) = sum(sum(It.*(x.*Ix+y.*Iy)));
B(2,1) = sum(sum(It.*(y.*Ix-x.*Iy)));
B(3,1) = sum(sum(It.*Ix));
B(4,1) = sum(sum(It.*Iy));
B = -B;

%%% Calculate the k parameters
C = inv(A)*B;
k = C';