
function res=bi_regressV2(X,Y,A,B,X1,Y1)
%
%use Friedman coding, should put the rotation matrix on left, change
%coordinate to a column. Also add the translation matrix just like time a
%matrix if we use augmented matrix which is used a lot in vizard coding
%e.g. [tx, ty] + [X, Y] = [1 0 tx; 0 1 ty; 0 0 1] * [X; Y; 1] = [X+tx; Y+ty; 1]
% X,Y,A,B,X1,Y1
% independent X, Y
% dependent A, B
% specific X1, Y1
%if we only have two pairs of X,Y,A,B, we get the perfect mapping just as
%we did using O and X1 in Mou & Zhang (2014), Angle=-h'




TemCov1= cov([A,X,B,Y]);

varA=TemCov1(1,1);
varX=TemCov1(2,2);
varB=TemCov1(3,3);
varY=TemCov1(4,4);
covAX=TemCov1(1,2);
covBY = TemCov1(3,4);
covAY = TemCov1(1,4);
covBX = TemCov1(3,2);


%X Y independent
beta1=(covAX + covBY)/(varX+varY); %X Y are independent variables
beta2=(covBX - covAY)/(varX+varY);

alpha1=mean(A)-beta1*mean(X)+beta2*mean(Y);
alpha2=mean(B)-beta2*mean(X)-beta1*mean(Y);

%Rp=[alpha1,alpha2]+

Scale=sqrt(beta1^2+beta2^2);
Angle=atan2(beta2,beta1)/pi()*180;

Ap=alpha1+beta1*X-beta2*Y;
Bp=alpha2+beta2*X+beta1*Y;

N=size(A,1);
unexplained=sum((A-Ap).^2+(B-Bp).^2);
explained=sum((Ap-mean(A)).^2+(Bp-mean(B)).^2);

Rsquared=explained/(explained+unexplained);
%Rsquared= 1-sum((A-Ap).^2+(B-Bp).^2)/((varA+varB)*(N-1));

[A, B, Ap, Bp];

DABsquared=unexplained;

Dmaxsquared=(N-1)*(varA+varB);

DIABsquared=DABsquared/Dmaxsquared; %which should be 1-Rsquared


% test it using P [1.378879998	0	1.157017697]


Ppx=alpha1+beta1*X1-beta2*Y1;
Ppy=alpha2+beta2*X1+beta1*Y1;



%output the regression and distortation variables
res.regression=[beta1,beta2, alpha1, alpha2,  Angle, Scale,Rsquared];
res.distortation=[DABsquared,Dmaxsquared, DIABsquared];
res.dependent=[Ppx,Ppy];
%close all
end