(* Compute the matrix of the projective transformation from    *)
(* the canonical quadrilateral to an arbitrary quadrilateral.  *)
fromCanonicaltoArbitrary[ q_ ] := 
Module [
{ M , A , P , b , l1 , l2 , l3 , res },
(
(* Compute the coefficients of the linear combination that  *)
(* defines the fourth vertex of the arbitrary quad in terms *)
(* of the other 3 points.                                   *)
 M = { 
       { q[[ 1 ]][[ 1 ]] , q[[ 2 ]][[ 1 ]] , q[[ 3 ]][[ 1 ]] } ,
       { q[[ 1 ]][[ 2 ]] , q[[ 2 ]][[ 2 ]] , q[[ 3 ]][[ 2 ]] } ,
       {                           1 ,                            1 ,                           1 }
     } ;

 b = { q[[ 4 ]][[ 1 ]] , q[[ 4 ]][[ 2 ]] ,  1 } ;

(* Test if the quadrangle (q1, q2, q3, q4) is convex    *)

{ l1 , l2 , l3 } = LinearSolve[ M , b ] ;
If[(l1 + l2 <= 1) || (l1<= 0) || (l2  ≤ 0), Print["** Warning, Nonconvex Quadrilateral **"],l1=l1]; 

(* Compute the matrix of the projective map that takes the *)
(* vertices of the canonical quad to the vertices of the   *)
(* arbitrary quad. The matrix is defined in terms of the   *)
(* basis defined by (1,0,1), (0,1,1), and (0,0,-1). So, we   *)
(* must perform a change of basis before multiplying this  *)
(* matrix by a point given in terms of the canonical basis *)
(* (1,0,0), (0,1,0), and (0,0,1).                          *) 
A = {
       { l1 q[[ 1 ]][[ 1 ]] , l2 q[[ 2 ]][[ 1 ]] , l3 q[[ 3 ]][[ 1 ]] } ,
       { l1 q[[ 1 ]][[ 2 ]] , l2 q[[ 2 ]][[ 2 ]] , l3 q[[ 3 ]][[ 2 ]] } ,
       {                              l1 ,                              l2 ,                               l3 } 
    } ;

(* Define the matrix that changes from the canonical basis *)
(* to the basis used to define the projective frame, which *)
(* is (1,0,1), (0,1,1), and (0,0,-1).                       *)
P = {
       {   1 ,   0 , 0  } ,
       {   0 ,   1 , 0  } ,
       { 1 ,1 , -1 } 
    } ;

res = A . P
) 
] ;

(* Show quadrilateral; remember, (1, 0), (0, 1), (0, 0), (1, 1) are mapped  *)
(*   to q1, q2, q3, q4                                                      *)

showquadrilateral[q_] :=
Module[
{qq1, qq2, qq3, qq4, res1, res2, res},
(qq1 = { q[[3]][[1]], q[[3]][[2]]};
qq2 = { q[[1]][[1]], q[[1]][[2]]};
qq3 = { q[[4]][[1]], q[[4]][[2]]};
qq4 = { q[[2]][[1]], q[[2]][[2]]};
res1 = Line[{qq1, qq2, qq3, qq4, qq1}];
res2 = {Graphics[Point[qq1]], Graphics[Point[qq2]], Graphics[Point[qq3]], Graphics[Point[qq4]]};
res = {res1, res2};
res)
];

(* Compute the matrix of the projective transformation from an  *)
(* arbitrary quadrilateral to another arbitrary quadrilateral.  *)
fromArbitrarytoArbitrary[ q1_ , q2_ ] := 
Module [
{ A1 , A2 , res },
(
A1 = fromCanonicaltoArbitrary[ q1 ] ;
A2 = fromCanonicaltoArbitrary[ q2 ] ;
res = A2 . Inverse[ A1 ] ;

res
)
] ;

(* Compute the image point under the projective transformation *)
(* that takes an arbitrary quadrilateral to another arbitrary  *)
(* quadrilateral.                                              *)
gFunc[ A_ , x_ , y_ ] := 
Module [
{  u , v , d , res },
(
u = A[[ 1 ]][[ 1 ]] x + A[[ 1 ]][[ 2 ]] y + A[[ 1 ]][[ 3 ]] ;
v = A[[ 2 ]][[ 1 ]] x + A[[ 2 ]][[ 2 ]] y + A[[ 2 ]][[ 3 ]] ;
d = A[[ 3 ]][[ 1 ]] x + A[[ 3 ]][[ 2 ]] y + A[[ 3 ]][[ 3 ]] ;
res = { u / d , v / d } ;

res
)
] ;

(*   computes and displays the image of a point from the canonical quadrilateral, in red   *)
(*  Also computes the line segment joining the x and y axes on which the map is undefined, in Orange *)
(*  and shows the original canonical  quadrangle can, in dotted lines   *)
(*  Shows the quadrangle (q1, q2, q3, q4) in green *)

showprojimage[q_,pt_] :=
Module[
{A, impt, t1, t2, t3, t4, t5, t6, xx, yy, p1, p2, can, res},
(
A = fromCanonicaltoArbitrary[q];
xx = pt[[1]]; yy = pt[[2]];
impt = gFunc[A, xx, yy];
p1 = {0, -A[[3]][[3]]/A[[3]][[2]]};
p2 = {-A[[3]][[3]]/A[[3]][[1]], 0};
(* Print["impt = ",impt]; *)
t1 = showquadrilateral[q];
t4 = Graphics[{Green, t1[[1]]}];
can =  {{1,0}, {0, 1}, {0, 0}, {1,1}};
t5 = showquadrilateral[can];
t6 = Graphics[{Dashed,t5[[1]]}];
t2 = Graphics[{PointSize[Medium],Red,Point[impt]}];
t3 = Graphics[{Orange,Line[{p1,p2}]}];
res = {t4,t1[[2]],t2,t3,t6};
res
)
];

ppp = {0,0.9};

(* A nonconvex polygon  *)

q = {{1, 0}, {0, 1}, {0, 0}, {0.4, 0.4}};

result = showprojimage[q,ppp];
Show[{result[[1]],result[[2]],result[[3]], result[[4]], result[[5]]}]

(*   computes and displays the image of a line segment from the canonical quadrilateral, in blue   *)
(*  Shows the original canonical  quadrangle can, in dashed lines   *)
(*  Shows the quadrangle (q1, q2, q3, q4) in green *)
(*  Also computes the line segment joining the x and y axes on which the map is undefined, in orange *)


showprojseg[q_,pt1_,pt2_] :=
Module[
{A, impt1, impt2, t1, t2, t3,t4, t5, t6, xt7, t8 , x1, yy1, xx2, yy2, p1, p2, can, t7, res},
(
A = fromCanonicaltoArbitrary[q];
xx1 = pt1[[1]]; yy1 = pt1[[2]];
impt1 = gFunc[A, xx1, yy1];
(* Print["impt = ",impt]; *)
xx2 = pt2[[1]]; yy2 = pt2[[2]];
impt2 = gFunc[A,xx2,yy2];
p1 = {0, -A[[3]][[3]]/A[[3]][[2]]};
p2 = {-A[[3]][[3]]/A[[3]][[1]], 0};
t1 = showquadrilateral[q];
t2 = Graphics[{PointSize[Medium],Red,Point[impt1]}];
t3 = Graphics[{PointSize[Medium],Blue, Point[impt2]}];
t4 = Graphics[{Blue, Line[{impt1,impt2}]}]; 
t5 = Graphics[{Orange,Line[{p1,p2}]}];
t6 = Graphics[{Green,t1[[1]]}];
can =  {{1,0}, {0, 1}, {0, 0}, {1,1}};
t7 = showquadrilateral[can];
t8 = Graphics[{Dashed, t7[[1]]}];
res = {t6,t1[[2]],t2, t3, t4, t5, t8};
res
)
];

pp1 = {0, 0.9}; pp2 = {1, 0.7};
result = showprojseg[q, pp1, pp2];
Show[{result[[1]], result[[2]], result[[3]], result[[4]], result[[5]], result[[6]], result[[7]]}]

(* A convex polygon *)
q = {{1, 0}, {0, 1}, {0, 0}, {2, 2}};

pp1 = {0, 0.2}; pp2 = {1, 0.7};
result = showprojseg[q, pp1, pp2];
Show[{result[[1]], result[[2]], result[[3]], result[[4]], result[[5]], result[[6]], result[[7]]}]

(* The h map - the double reflection *)
hFunc[ x_ , y_ ] :=
Module [
{ res  } ,
(
res = { 1 - x , -y } ;
res
)
] ;

(* A rotation map - a general one for all gluing domains *)
rFunc[ x_ , y_ , t_ ] :=
Module [
{ cost , sint , u , v , res } ,
(
cost = Cos[ t ] ;
sint = Sin[ t ] ;
u = x cost - y sint ;
v = x sint + y cost ;
res = { u , v } ;
res
)
] ;

(* Quadrilateral 1 - corresponding to the p-domain Omega_u  *)
qu = { { Cos[ -2Pi / 5 ] ,Sin[ -2Pi / 5 ] } , { Cos[ 2Pi / 5 ] ,Sin[ 2Pi / 5 ] } , { 0 , 0 } ,{ 1 , 0 } } ;

(* Quadrilateral 2 - corresponding to the p-domain Omega_v  *)
qv = { { Cos[ -2Pi / 6 ] ,Sin[ -2Pi / 6 ] } , { Cos[ 2Pi / 6 ] ,Sin[ 2Pi / 6 ] } , { 0 , 0 } ,{ 1 , 0 }  } ;

(* Quadrilateral 3 - corresponding to the p-domain Omega_w  *)
qw = { { Cos[ -2Pi / 7 ] ,Sin[ -2Pi / 7 ] } , { Cos[ 2Pi / 7 ] ,Sin[ 2Pi / 7 ] } , { 0 , 0 } ,{ 1 , 0 } } ;

(* Canonical quadrilateral *)
qq = { { Cos[ -2Pi / 6 ] ,Sin[ -2Pi / 6 ] } , { Cos[ 2Pi / 6 ] ,Sin[ 2Pi / 6 ] } , { 0 , 0 } ,{ 1 , 0 }  } ;

(* Compute the projective transformation matrix from qu to qq *)
Mu = fromArbitrarytoArbitrary[ qu , qq ] ;

(* Compute the projective transformation matrix from qq to qu *)
MuInv = fromArbitrarytoArbitrary[ qq , qu ] ;

(* Compute the projective transformation matrix from qv to qq *)
Mv = fromArbitrarytoArbitrary[ qv , qq ] ;

(* Compute the projective transformation matrix from qq to qv *)
MvInv = fromArbitrarytoArbitrary[ qq , qv ] ;

(* Compute the projective transformation matrix from qw to qq *)
Mw = fromArbitrarytoArbitrary[ qw , qq ] ;

(* Compute the projective transformation matrix from qq to qw *)
MwInv = fromArbitrarytoArbitrary[ qq , qw ] ;

(* Compute the projective transformation from qu to qq *)
fromUtoQ[ x_ , y_ ] := gFunc[ Mu , x , y ] ;

(* Compute the projective transformation from qq to qu *)
fromQtoU[ x_ , y_ ] := gFunc[ MuInv , x , y ] ;

(* Compute the projective transformation from qv to qq *)
fromVtoQ[ x_ , y_ ] := gFunc[ Mv , x , y ] ;

(* Compute the projective transformation from qq to qv *)
fromQtoV[ x_ , y_ ] := gFunc[ MvInv , x , y ] ;

(* Compute the projective transformation from qw to qq *)
fromWtoQ[ x_ , y_ ] := gFunc[ Mw , x , y ] ;

(* Compute the projective transformation from qq to qw *)
fromQtoW[ x_ , y_ ] := gFunc[ MwInv , x , y ] ;

(* Testing the transformations *)
N[ fromUtoQ[ 0 , 0 ] ]

N[ fromUtoQ[ 1 , 0 ] ]

N[ fromUtoQ[ Cos[ -2Pi / 5 ], Sin[ -2Pi / 5 ] ] ]

N[ fromUtoQ[ Cos[ 2Pi / 5 ], Sin[ 2Pi / 5 ] ] ]

N[ fromQtoU[ 0 , 0 ] ]

N[ fromQtoU[ 1 , 0 ] ]

N[ fromQtoU[ Cos[ -Pi / 3 ], Sin[ -Pi / 3 ] ] ]

N[ fromQtoU[ Cos[ Pi / 3 ], Sin[ Pi / 3 ] ] ]

N[ fromVtoQ[ 0 , 0 ] ]

N[ fromVtoQ[ 1 , 0 ] ]

N[ fromVtoQ[ Cos[ -Pi / 3 ], Sin[ -Pi / 3 ] ] ]

N[ fromVtoQ[ Cos[ Pi / 3 ], Sin[ Pi / 3 ] ] ]

N[ fromQtoV[ 0 , 0 ] ]

N[ fromQtoV[ 1 , 0 ] ]

N[ fromQtoV[ Cos[-Pi/3  ] , Sin[ -Pi / 3 ] ] ]

N[ fromQtoV[ Cos[ Pi/3  ] , Sin[ Pi / 3 ] ] ]

N[ fromWtoQ[ 0 , 0 ] ]

N[ fromWtoQ[ 1 , 0 ] ]

N[ fromWtoQ[ Cos[-2Pi/7  ] , Sin[ -2Pi / 7 ] ] ]

N[ fromWtoQ[ Cos[ 2Pi/7  ] , Sin[ 2Pi / 7 ] ] ]

N[ fromQtoW[ 0 , 0 ] ]

N[ fromQtoW[ 1 , 0 ] ]

N[ fromQtoW[ Cos[-Pi/3  ] , Sin[ -Pi / 3 ] ] ]

N[ fromQtoW[ Cos[ Pi/3  ] , Sin[ Pi / 3 ] ] ]

(* Testing the cocycle condition - Assume that w precedes v in *)
(* a counterclockwise traversal of lk( u , K ).                *)

(* Transition function from Omega_{uv} to Omega_{vu} *)
phiVU[ x_ , y_ ] :=
Module[
{ r , s , res } ,
(
{ r , s } = rFunc[ x , y , -2Pi / 5 ] ;
{ r , s } = fromUtoQ[ r , s ] ;
{ r , s } = hFunc[ r , s ] ;
{ r , s } = fromQtoV[ r , s ] ;
res = { r , s } ;
res
)
] ;

(* Transition function from Omega_{uw} to Omega_{wu} *)
phiWU[ x_ , y_ ] :=
Module[
{ r , s , res } ,
(
{ r , s } = fromUtoQ[ x , y ] ;
{ r , s } = hFunc[ r , s ] ;
{ r , s } = fromQtoW[ r , s ] ;
{ r , s } = rFunc[ r , s ,2Pi / 7 ] ;
res = { r , s } ;
res
)
] ;

(* Transition function from Omega_{vu} to Omega_{uv} *)
phiUV[ x_ , y_ ] :=
Module[
{ r , s , res } ,
(
{ r , s } = fromVtoQ[ x , y ] ;
{ r , s } = hFunc[ r , s ] ;
{ r , s } = fromQtoU[ r , s ] ;
{ r , s } = rFunc[ r , s , 2Pi / 5 ] ;
res = { r , s } ;
res
)
] ;

(* Transition function from Omega_{vw} to Omega_{wv} *)
phiWV[ x_ , y_ ] :=
Module[
{ r , s , res } ,
(
{ r , s } = rFunc[ x , y , -Pi / 3 ] ;
{ r , s } = fromVtoQ[ r , s ] ;
{ r , s } = hFunc[ r , s ] ;
{ r , s } = fromQtoW[ r , s ] ;
res = { r , s } ;
res
)
] ;

(* Transition function from Omega_{wu} to Omega_{uw} *)
phiUW[ x_ , y_ ] :=
Module[
{ r , s , res } ,
(
{ r , s } = rFunc[ x , y , -2Pi / 7 ] ;
{ r , s } = fromWtoQ[ r , s ] ;
{ r , s } = hFunc[ r , s ] ;
{ r , s } = fromQtoU[ r , s ] ;
res = { r , s } ;
res
)
] ;

(* Transition function from Omega_{wv} to Omega_{vw} *)
phiVW[ x_ , y_ ] :=
Module[
{ r , s , res } ,
(
{ r , s } = fromWtoQ[ x , y ] ;
{ r , s } = hFunc[ r , s ] ;
{ r , s } = fromQtoV[ r , s ] ;
{ r , s } = rFunc[ r , s , Pi / 3 ] ;
res = { r , s } ;
res
)
] ;

(* Testing the transition functions *)

N[ phiVU[ 0,0] ]

N[ phiVU[ 1,0] ]

N[ phiVU[ Cos[2Pi/5],Sin[2Pi/5]] ]

N[ phiWU[ 0,0] ]

N[ phiWU[ 1,0] ]

N[ phiWU[ Cos[2Pi/5],Sin[2Pi/5]] ]

N[ phiUV[ 0,0] ]

N[ phiUV[ 1,0] ]

N[ phiUV[ Cos[ Pi / 3 ],Sin[ Pi / 3 ] ] ]

N[ phiWV[ 0,0] ]

N[ phiWV[ 1,0] ]

N[ phiWV[ Cos[ Pi / 3 ],Sin[ Pi / 3 ]] ]

N[ phiUW[ 0,0] ]

N[ phiUW[ 1,0] ]

N[ phiUW[ Cos[2Pi/7],Sin[2Pi/7] ] ]

N[ phiVW[ 0,0] ]

N[ phiVW[ 1,0] ]

N[ phiVW[ Cos[2Pi/7],Sin[2Pi/7] ] ]

(* Show that the projective transformation does not *)
(* satisfy the cocycle condition.                   *)

{ x1 , y1 } = phiVU[ 0.5 , 0.5 ]

{ x2 , y2 } = phiWU[ 0.5 , 0.5 ]

{ x3 , y3 } = phiVW[ x2 , y2 ]

N[ Mu ]

N[ MuInv ]

N[ Mv ]

N[ MvInv ]

N[ Mw ]

N[ MwInv ]

N[ Mu ] . N[ MuInv ]

N[ Mw ] . N[ MwInv ]