% Define transformation
M = [ 2 3 1 ; 3 1 0; 1 2 1]
M =[  0.7071     0          1.0 ;
         0.7071    1.4142     0 ;
         1.0000    2.0        1.0]   
   
 % Define xy coordinates of unit square
 xy =[
     0     1     1     0;
     0     0     1     1]
 xy_hom = [xy ; ones(1,size(xy,2))]
 
 % Calculate transformed coordinates of unit square
 xynew_hom = M*xy_hom
 xynew = xynew_hom(1:2,:) ./ repMat(xynew_hom(3,:),[2,1])
 
 % plot orig (red) and transformed (blue)
 figure;
 plot([xy(1,:),xy(1,1)],[xy(2,:),xy(2,1)],'r','Linewidth',2);
 hold on
 plot([xynew(1,:),xynew(1,1)],[xynew(2,:),xynew(2,1)],'b','Linewidth',2);
 axis equal
 
 % Analyze transform
 M_nonproj = M; M_nonproj(3,1:2)=0 % define non-projective part of M
 M_proj = M*inv(M_nonproj)  %we want M_proj*M_nonproj = M
 (M_proj*M_nonproj) - M  % check
 
 M_affine = M_nonproj; M_affine(1:2,3) = 0 % define affine part of M
 M_tran = M_nonproj*inv(M_affine)  % we want M_tran*M_affine = M_nonproj
 (M_proj*M_tran*M_affine) - M  % check
 
 % Guess Rotation 45 deg
 M_rot = [sqrt(2)/2 -sqrt(2)/2  0; sqrt(2)/2 sqrt(2)/2 0 ; 0 0 1]
 M_base = inv(M_rot)*M_affine % M_rot*M_base = M_affine
 
 % Final check
M_proj*M_tran*M_rot*M_base - M
 