<p>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.</p>

<p>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).</p>

<p>Let <span class="math-container">$M$</span> be the transform matrix, then</p>

<p><span class="math-container">$$
\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}
$$</span>
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 <span class="math-container">$(X_4,Y_4)$</span>:</p>

<p><span class="math-container">$$
\begin{align}
H_4 &amp;= M (1,1,1) \\
\frac{H_4}{H_4.z} &amp;= (X_4,Y_4,1)
\end{align}
$$</span>
Note that I have to divide the homogeneous part for the transformed point.</p>

<p>Solve this linear system you can get an expression for your matrix <span class="math-container">$M$</span>. Perform <span class="math-container">$M (f_x,f_y,1)$</span> you will see how it transforms the point. I believe this solves your Q2, Q3, Q5 and Q7.</p>

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

<hr>

<p>Appendix: maple code to verify the answer</p>

<pre><code>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
</code></pre>

CC BY-SA 4.0
2019-02-20T21:55:41.273Z
2019-02-20T21:55:41.273Z
TheBusyTypist
<p>I've spent about 2 days trying to understand this piece of code (from <a href="https://github.com/nilx/megawave/blob/master/src/image/domain/fproj.c" rel="nofollow noreferrer">here</a>) that applies an affine/projective transformation to an image. I will link bits of the code that I'm struggling to understand and add comments on what I think it's trying to do. In particular, the author's original comments I'll label as <code>//</code>. </p>

<pre><code>//      COMPUTE NEW BASIS
// X1, Y1 : upleft corner
// X2, Y2 : upright corner
// X3, Y3 : downleft corner 
// sx     : x-size of output image
// sy     : y-size of output image
x12 = (X2-X1)/(float)(*sx); 
y12 = (Y2-Y1)/(float)(*sx);
x13 = (X3-X1)/(float)(*sy);
y13 = (Y3-Y1)/(float)(*sy);
</code></pre>

<p>Q1. I can see that this part as the comment says, finds a new normalized basis but I'm unsure exactly why we need to find one to model an affine transformation?</p>

<pre><code>// x4, y4 : downright corner (for projective transform)

if (y4) 
{ 
    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));
    a = (yy-1.0)/(1.0-xx-yy);
    b = (xx-1.0)/(1.0-xx-yy);
} 
else 
{
    a=b=0.0;
}
</code></pre>

<p>After doing some internet search, I manage to find equations from <a href="http://paulbourke.net/geometry/pointlineplane/" rel="nofollow noreferrer">here</a> that the <code>xx</code> and <code>yy</code> equations solve for the "intersection point of two line segments in 2 dimensions". An image will be linked below.</p>

<p><strong>START IMAGE.</strong>
<a href="https://i.sstatic.net/ANCv2m.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/ANCv2m.png" alt="Intersection of two lines"></a></p>

<p><strong>END IMAGE.</strong></p>

<p>Q2. I don't understand how he uses the <code>xx</code> and <code>yy</code> equations to solve for a &amp; b? </p>

<p>Q3. I am unsure of what <code>a</code> and <code>b</code> are. My best guess is that they are the bottom row of a projective transformation matrix shown below:</p>

<p>$$
\begin{bmatrix}
a_{11} &amp; a_{12} &amp; a_{13} \\
a_{21} &amp; a_{22} &amp; a_{23} \\
a &amp; b &amp; 1 \\
\end{bmatrix}
$$</p>

<p>This is because the code sets <code>a = b = 0.0</code> if we're not doing a projective transform, which should then make it an affine transformation matrix.</p>

<p>Now is where the main loop happens:</p>

<pre><code>for (x=0;x&lt;sx;x++) 
    for (y=0;y&lt;sy;y++) {
        fx = (float)x + 0.5;
        fy = (float)y + 0.5;
        d = a*fx/(float)(*sx)+b*fy/(float)(*sy)+1.0;
        xx = (a+1.0)*fx/d;
        yy = (b+1.0)*fy/d;
        xp = X1 + xx*x12 + yy*x13;
        yp = Y1 + xx*y12 + yy*y13;
    }
</code></pre>

<p>Q4. Why does he add 0.5 to the original x and y coordinate? My best guess is that the code's trying to map the centre of each pixel rather than the edge to the new image:</p>

<p><a href="https://i.sstatic.net/vnPXU.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/vnPXU.png" alt="Pixels"></a></p>

<p>Q5. In </p>

<pre><code>d = a*fx/(float)(*sx)+b*fy/(float)(*sy)+1.0;
xx = (a+1.0)*fx/d;
yy = (b+1.0)*fy/d;
</code></pre>

<p>what are they trying to accomplish here? My best guess is that <code>d</code> is the third component of the output homogeneous coordinate:</p>

<p>$$
\begin{bmatrix}
xx\\
yy\\
d\\
\end{bmatrix}
$$</p>

<p>calculated from the dot product of the third row of the (now normalised) transformation matrix from earlier and the input homogeneous vector:</p>

<p>$$
\begin{bmatrix}
\frac{a}{sx} &amp; \frac{b}{sy} &amp; 1 \\    
\end{bmatrix}
\begin{bmatrix}
fx\\
fy\\
1\\
\end{bmatrix}
=
a*fx/sx+b*fy/sy+1.0
$$</p>

<p>Then the division</p>

<pre><code>xx = (a+1.0)*fx/d;
yy = (b+1.0)*fy/d;
</code></pre>

<p>would just transform the homogeneous coordinates back to Euclidean coordinates, corresponding to the output image pixel location. I noted that if <code>a=b=0.0</code> (so not projective), then <code>(a+1.0)=(b+1.0)=d=1.0</code>. I'm however not sure about why we have to multiply by <code>(a+1.0)</code> or <code>(b+1.0)</code>. But again, I do not know if I'm on the right track or not as I do not have anyone to help me with this code.</p>

<p>Q7. The final operation</p>

<pre><code>xp = X1 + xx*x12 + yy*x13;
yp = Y1 + xx*y12 + yy*y13;
</code></pre>

<p>I suppose is literally the dot product of the first two rows of the transformation matrix with the Euclidean coordinates as so:</p>

<p>$$
\begin{bmatrix}
x12 &amp; x13 &amp; X_{1} \\ 
y12 &amp; y13 &amp; Y_{1}  
\end{bmatrix}
\begin{bmatrix}
xx\\
yy\\
1\\
\end{bmatrix}
=
\begin{bmatrix}
xp\\
yp\\
\end{bmatrix}
$$</p>

<p>Which brings me to my question, how does this achieve an affine transformation? </p>

<p>I suspect my lack of formal training could be why I'm not understanding how the code models an affine transformation. I'd appreciate some light shed on this piece of code!</p>

<p>Cheers!</p>

CC BY-SA 4.0
2018-09-11T08:23:43.490Z
An affine transformation is is any transformation that preserves points, lines and planes. This can be translation, skewing and rotation - and combinations thereof.
beyond
CC BY-SA 4.0
2018-09-11T09:00:40.970Z
Hi @beyond, thanks for your input. Are you able to shed some light on the code?
theunseen
CC BY-SA 4.0
2018-09-11T12:57:41.527Z
To be honest, the code seems a bit elaborate for something as simple as applying an affine transformation to an image. The projection bit I don't get, so an example of input/expected output/actual output would be ideal.
beyond
<transformations><algorithm><affine-transformations><linear-algebra>
Deciphering Affine/Projective Transformation Code
theunseen
# Deciphering Affine/Projective Transformation Code

(Asked by: username_0 on 2018-09-10)

I've spent about 2 days trying to understand this piece of code (from [here](https://github.com/nilx/megawave/blob/master/src/image/domain/fproj.c)) that applies an affine/projective transformation to an image. I will link bits of the code that I'm struggling to understand and add comments on what I think it's trying to do. In particular, the author's original comments I'll label as `//`.

```
//      COMPUTE NEW BASIS
// X1, Y1 : upleft corner
// X2, Y2 : upright corner
// X3, Y3 : downleft corner 
// sx     : x-size of output image
// sy     : y-size of output image
x12 = (X2-X1)/(float)(*sx); 
y12 = (Y2-Y1)/(float)(*sx);
x13 = (X3-X1)/(float)(*sy);
y13 = (Y3-Y1)/(float)(*sy);

```

Q1. I can see that this part as the comment says, finds a new normalized basis but I'm unsure exactly why we need to find one to model an affine transformation?

```
// x4, y4 : downright corner (for projective transform)

if (y4) 
{ 
    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));
    a = (yy-1.0)/(1.0-xx-yy);
    b = (xx-1.0)/(1.0-xx-yy);
} 
else 
{
    a=b=0.0;
}

```

After doing some internet search, I manage to find equations from [here](http://paulbourke.net/geometry/pointlineplane/) that the `xx` and `yy` equations solve for the "intersection point of two line segments in 2 dimensions". An image will be linked below.

**START IMAGE.**
[![Intersection of two lines](https://i.sstatic.net/ANCv2m.png)](https://i.sstatic.net/ANCv2m.png)

**END IMAGE.**

Q2. I don't understand how he uses the `xx` and `yy` equations to solve for a & b?

Q3. I am unsure of what `a` and `b` are. My best guess is that they are the bottom row of a projective transformation matrix shown below:

$$
\begin{bmatrix}
a\_{11} & a\_{12} & a\_{13} \\
a\_{21} & a\_{22} & a\_{23} \\
a & b & 1 \\
\end{bmatrix}
$$

This is because the code sets `a = b = 0.0` if we're not doing a projective transform, which should then make it an affine transformation matrix.

Now is where the main loop happens:

```
for (x=0;x<sx;x++) 
    for (y=0;y<sy;y++) {
        fx = (float)x + 0.5;
        fy = (float)y + 0.5;
        d = a*fx/(float)(*sx)+b*fy/(float)(*sy)+1.0;
        xx = (a+1.0)*fx/d;
        yy = (b+1.0)*fy/d;
        xp = X1 + xx*x12 + yy*x13;
        yp = Y1 + xx*y12 + yy*y13;
    }

```

Q4. Why does he add 0.5 to the original x and y coordinate? My best guess is that the code's trying to map the centre of each pixel rather than the edge to the new image:

[![Pixels](https://i.sstatic.net/vnPXU.png)](https://i.sstatic.net/vnPXU.png)

Q5. In

```
d = a*fx/(float)(*sx)+b*fy/(float)(*sy)+1.0;
xx = (a+1.0)*fx/d;
yy = (b+1.0)*fy/d;

```

what are they trying to accomplish here? My best guess is that `d` is the third component of the output homogeneous coordinate:

$$
\begin{bmatrix}
xx\\
yy\\
d\\
\end{bmatrix}
$$

calculated from the dot product of the third row of the (now normalised) transformation matrix from earlier and the input homogeneous vector:

$$
\begin{bmatrix}
\frac{a}{sx} & \frac{b}{sy} & 1 \\
\end{bmatrix}
\begin{bmatrix}
fx\\
fy\\
1\\
\end{bmatrix}
=
a\*fx/sx+b\*fy/sy+1.0
$$

Then the division

```
xx = (a+1.0)*fx/d;
yy = (b+1.0)*fy/d;

```

would just transform the homogeneous coordinates back to Euclidean coordinates, corresponding to the output image pixel location. I noted that if `a=b=0.0` (so not projective), then `(a+1.0)=(b+1.0)=d=1.0`. I'm however not sure about why we have to multiply by `(a+1.0)` or `(b+1.0)`. But again, I do not know if I'm on the right track or not as I do not have anyone to help me with this code.

Q7. The final operation

```
xp = X1 + xx*x12 + yy*x13;
yp = Y1 + xx*y12 + yy*y13;

```

I suppose is literally the dot product of the first two rows of the transformation matrix with the Euclidean coordinates as so:

$$
\begin{bmatrix}
x12 & x13 & X\_{1} \\
y12 & y13 & Y\_{1}
\end{bmatrix}
\begin{bmatrix}
xx\\
yy\\
1\\
\end{bmatrix}
=
\begin{bmatrix}
xp\\
yp\\
\end{bmatrix}
$$

Which brings me to my question, how does this achieve an affine transformation?

I suspect my lack of formal training could be why I'm not understanding how the code models an affine transformation. I'd appreciate some light shed on this piece of code!

Cheers!

**Tags:** transformations, algorithm, affine-transformations, linear-algebra

> **Comments:**
> **username_1**: An affine transformation is is any transformation that preserves points, lines and planes. This can be translation, skewing and rotation - and combinations thereof.  
> **username_0**: Hi @username_1, thanks for your input. Are you able to shed some light on the code?  
> **username_1**: To be honest, the code seems a bit elaborate for something as simple as applying an affine transformation to an image. The projection bit I don't get, so an example of input/expected output/actual output would be ideal.  

---

## 1 answers

### Answer by username_2 on 2019-02-20. Score: 2 | Accepted answer

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) &= (X\_1,Y\_1,1) \\
M (1,0,1) &= (X\_2,Y\_2,1) \\
M (0,1,1) &= (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 &= M (1,1,1) \\
\frac{H\_4}{H\_4.z} &= (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:=<m11,m21,0|m12,m22,0|m13,m23,1>;
# Build equations
E:=[convert(m.<0,0,1> - <X1,Y1,1>,list),convert(m.<1,0,1> - <X2,Y2,1>,list), convert(m.<0,1,1> - <X3,Y3,1>,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:=<m11,m21,a|m12,m22,b|m13,m23,1>;
# Compute transformed points
U0:=p.<0,0,1>;
U1:=p.<1,0,1>;
U2:=p.<1,1,1>;
U2:=p.<0,1,1>;
U3:=p.<1,1,1>;
# Build equations
E:=[convert(U0/U0[3]-<X1,Y1,1>,list),convert(U1/U1[3]-<X2,Y2,1>,list),convert(U2/U2[3]-<X3,Y3,1>,list),convert(U3/U3[3]-<X4,Y4,1>,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

```


---

