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'



figure
plot (X, Y,'k')
text(X(1)+0.1,Y(1),['O',char(39)],'Color', 'k');%['O', char(39)] means O'
hold on
plot (A, B,'r')

text(A(1)+0.1,B(1),'O', 'Color', 'r');%['O', char(39)] means O'
hold on
%now let's use Friedman equations to predict M (or X Y) using R (or A B)
%TemCov1= cov([A(1:2),X(1:2),B(1:2),Y(1:2)])
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

plot (Ap, Bp,'b')
hold on
% test it using P [1.378879998	0	1.157017697]


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


plot (X1, Y1,'k')
text(X1+0.1,Y1,'p','Color', 'k');%['O', char(39)] means O'
hold on
plot (Ppx, Ppy,'b')
text(Ppx+0.1,Ppy,['p',char(39)],'Color', 'b');%['O', char(39)] means O'
%hold on

grid on
pbaspect([1 1 1])
axislim=20;
xlim([-axislim axislim])
ylim([-axislim axislim])

%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