function [px,py]=transform_point(Ax,Ay,Bx,By,Px,Py,ax,ay,bx,by)

% [Ax,Ay] and [Bx,By] are the initial coordinate system reference
% [ax,ay] and [bx,by] are the new coordinate system reference

% [Px,Py] is the point to be transformed
% [px,py] is the transformed coordinates of [Px,Py]

% Compute coordinate system scale factor

F = sqrt((bx-ax)^2+(by-ay)^2)/sqrt((Bx-Ax)^2+(By-Ay)^2);

% Compute angles that define coordinate system rotation

Theta=atan((By-Ay)/(Bx-Ax));
theta=atan((by-ay)/(bx-ax));

% Correct returned angle taking into account correct quadrant

if (((Theta > 0) && ((By-Ay) < 0)) || ((Theta < 0) && ((By-Ay) > 0)))
    Theta=Theta+pi;
end

if (((theta > 0) && ((by-ay) < 0)) || ((theta < 0) && ((by-ay) > 0)))
    theta=theta+pi;
end

px = F*((Px-Ax)*cos(theta-Theta)-(Py-Ay)*sin(theta-Theta))+ax;
py = F*((Px-Ax)*sin(theta-Theta)+(Py-Ay)*cos(theta-Theta))+ay;

figure;
hold on;
plot(Ax,Ay,'*');
plot(Bx,By,'o');
plot(Px,Py,'+');
title('Initial coordinate system');
legend('A','B','P')

figure;
hold on;
plot(ax,ay,'*');
plot(bx,by,'o');
plot(px,py,'+');
title('New coordinate system');
legend('a','b','p')

end