Deciphering Affine/Projective Transformation Code
Your understanding of the matrix structure in Q3 is correct. This code just does not construct a matrix explicitly and the matrix multiplication is applied implicitly. I think this part might cause your confusion.
Instead of deciphering the code, I would rather derive the transform and compare it with the code.

The affine (6 degrees of freedom) and projective matrix (8 dof) map points to points. So you can build a system of equations for each pair of points in different frames. You will need to solve this system for the components of the matrix (Q1).

Let $M$ be the transform matrix, then

$$
\begin{align}
M (0,0,1) &amp;= (X_1,Y_1,1) \\
M (1,0,1) &amp;= (X_2,Y_2,1) \\
M (0,1,1) &amp;= (X_3,Y_3,1) \\
\end{align}
$$
Each of above lines gives you two equations. Now you have 6 equations and the affine matrix can be derived.
For the projective matrix, we need extra two equations using the point $(X_4,Y_4)$:

$$
\begin{align}
H_4 &amp;= M (1,1,1) \\
\frac{H_4}{H_4.z} &amp;= (X_4,Y_4,1)
\end{align}
$$
Note that I have to divide the homogeneous part for the transformed point.

Solve this linear system you can get an expression for your matrix $M$. Perform $M (f_x,f_y,1)$ you will see how it transforms the point. I believe this solves your Q2, Q3, Q5 and Q7.

For Q4, I think it is just a convention of pixel position and your understanding should be correct.



Appendix: maple code to verify the answer

with(ListTools):
# affine matrix
m:=&lt;m11,m21,0|m12,m22,0|m13,m23,1&gt;;
# Build equations
E:=[convert(m.&lt;0,0,1&gt; - &lt;X1,Y1,1&gt;,list),convert(m.&lt;1,0,1&gt; - &lt;X2,Y2,1&gt;,list), convert(m.&lt;0,1,1&gt; - &lt;X3,Y3,1&gt;,list)];
# E=[[m13 - X1, m23 - Y1, 0], [m11 + m13 - X2, m21 + m23 - Y2, 0],[m12 + m13 - X3, m22 + m23 - Y3, 0]]
# Solve E=0
sol:=eliminate(Flatten(E),indets(m));
# Plug the solution back to m:
subs(sol[1],m);
                           [-X1 + X2    -X1 + X3    X1]
                           [                          ]
                           [-Y1 + Y2    -Y1 + Y3    Y1]
                           [                          ]
                           [   0           0        1 ]


# projective matrix
p:=&lt;m11,m21,a|m12,m22,b|m13,m23,1&gt;;
# Compute transformed points
U0:=p.&lt;0,0,1&gt;;
U1:=p.&lt;1,0,1&gt;;
U2:=p.&lt;1,1,1&gt;;
U2:=p.&lt;0,1,1&gt;;
U3:=p.&lt;1,1,1&gt;;
# Build equations
E:=[convert(U0/U0[3]-&lt;X1,Y1,1&gt;,list),convert(U1/U1[3]-&lt;X2,Y2,1&gt;,list),convert(U2/U2[3]-&lt;X3,Y3,1&gt;,list),convert(U3/U3[3]-&lt;X4,Y4,1&gt;,list)];
# Sove E=0
solp:=eliminate(Flatten(E),indets(p));
# copy their code for xx and yy:
xx:=((X4-X1)*(Y3-Y1)-(Y4-Y1)*(X3-X1))/((X2-X1)*(Y3-Y1)-(Y2-Y1)*(X3-X1));
yy:=((X4-X1)*(Y2-Y1)-(Y4-Y1)*(X2-X1))/((X3-X1)*(Y2-Y1)-(Y3-Y1)*(X2-X1));
# Verify our a and b, compare ours with their expressing using xx and yy
subs(solp[1],a);
subs(solp[1],b);
simplify(subs(solp[1],a)-(yy-1)/(1-xx-yy));
# It is 0 as expected
simplify(subs(solp[1],b)-(xx-1)/(1-xx-yy));
# It is also 0