%!

% a matrix is an array [[...][...][...]] of columns

%v M -> v.M
% [v0 v1 v2][M0 M1 M2]
/vector-matrix-mul { 1 dict begin
/M exch def
/v exch def
[
  0 1 2 { /i exch def
    v 0 get M i get 0 get mul
    v 1 get M i get 1 get mul add
    v 2 get M i get 2 get mul add
  } for
]
end } def


% M N -> M.N
% [A0 A1 A2][B0 B1 B2]
/matrix-matrix-mul { 1 dict begin
/B exch def
/A exch def
[
  0 1 2 { /j exch def
  [
    0 1 2 { /i exch def
      A 0 get i get B j get 0 get mul
      A 1 get i get B j get 1 get mul add
      A 2 get i get B j get 2 get mul add
    } for
  ]
  } for
]


end } def

% T -> det(T)
/3x3determinant { 1 dict begin
/T exch def
  T 0 get 0 get 
  T 1 get 1 get mul
  T 2 get 2 get mul

  T 0 get 1 get 
  T 1 get 2 get mul
  T 2 get 0 get mul add

  T 0 get 2 get 
  T 1 get 0 get mul
  T 2 get 1 get mul add

  T 0 get 2 get 
  T 1 get 1 get mul
  T 2 get 0 get mul sub

  T 0 get 1 get 
  T 1 get 0 get mul
  T 2 get 2 get mul sub

  T 0 get 0 get 
  T 1 get 2 get mul
  T 2 get 1 get mul sub
end } def

% T -> T^{-1}
/3x3inverse { 1 dict begin
/T exch def
/D T 3x3determinant def
[
  [
    T 1 get 1 get T 2 get 2 get mul
    T 1 get 2 get T 2 get 1 get mul sub
	D div

    T 0 get 1 get T 2 get 2 get mul
    T 0 get 2 get T 2 get 1 get mul sub
	D div neg

    T 0 get 1 get T 1 get 2 get mul
    T 0 get 2 get T 1 get 1 get mul sub
	D div
  ]
  [
    T 1 get 0 get T 2 get 2 get mul
    T 1 get 2 get T 2 get 0 get mul sub
	D div neg

    T 0 get 0 get T 2 get 2 get mul
    T 0 get 2 get T 2 get 0 get mul sub
	D div 

    T 0 get 0 get T 1 get 2 get mul
    T 0 get 2 get T 1 get 0 get mul sub
	D div neg
  ]
  [
    T 1 get 0 get T 2 get 1 get mul
    T 1 get 1 get T 2 get 0 get mul sub
	D div

    T 0 get 0 get T 2 get 1 get mul
    T 0 get 1 get T 2 get 0 get mul sub
	D div neg

    T 0 get 0 get T 1 get 1 get mul
    T 0 get 1 get T 1 get 0 get mul sub
	D div
  ]
]
end } def

% x y T - > transform in 2D
/projective-transform { 1 dict begin
/T exch def
/y exch def
/x exch def
/x' 
  T 0 get 0 get x mul 
  T 0 get 1 get y mul add
  T 0 get 2 get add 
def
/y' 
  T 1 get 0 get x mul 
  T 1 get 1 get y mul add
  T 1 get 2 get add 
def
/z' 
  T 2 get 0 get x mul 
  T 2 get 1 get y mul add
  T 2 get 2 get add 
def
x' z' div
y' z' div 
end } def

% [P Q R S] -> 3x3 projective matrix
/projective-matrix { 1 dict begin
aload pop % P Q R S
aload pop 
/ys exch def
/xs exch def

aload pop 
/yr exch def
/xr exch def

aload pop 
/yq exch def
/xq exch def

aload pop 
/yp exch def
/xp exch def

% find cp cq cr  
/C  [0 1 1]   [[0 1 1][0 0 1][1 1 1]] 3x3inverse vector-matrix-mul def
/C' [xs ys 1] [[xp xq xr][yp yq yr][1 1 1]] 3x3inverse vector-matrix-mul def
/M  [[C 0 get 0 mul 
      C 1 get 1 mul 
      C 2 get 1 mul]
     [C 0 get 0 mul 
      C 1 get 0 mul 
      C 2 get 1 mul]
     [C 0 get 
      C 1 get 
      C 2 get]
    ] def
/M'  [[C' 0 get xp mul 
      C' 1 get xq mul 
      C' 2 get xr mul]
     [C' 0 get yp mul 
      C' 1 get yq mul 
      C' 2 get yr mul]
     [C' 0 get 
      C' 1 get 
      C' 2 get]
    ] def

M 3x3inverse M' matrix-matrix-mul 

end } def

% ---------------------------------------------------------------

/S 72 3 mul def
S dup scale
1 S div setlinewidth

1 1 translate

/T [[0 0][1 0][1 1][0 1]] projective-matrix def

/N 10 def

% square

% transformed

gsave
1 0 0 setrgbcolor
newpath
0 1 N { /i exch def
/I i N div def
0 I T projective-transform moveto
1 I T projective-transform lineto
I 0 T projective-transform moveto
I 1 T projective-transform lineto
} for
closepath
stroke
grestore

/T [[0 0][1 0][0.75 0.5][0.25 0.5]] projective-matrix def

newpath
0 1 N { /i exch def
/I i N div def
0 I T projective-transform moveto
1 I T projective-transform lineto
I 0 T projective-transform moveto
I 1 T projective-transform lineto
} for
closepath
stroke