How to find the (projective) inscribed circle in a quadrilateral

This post answers a question raised in Extraordinary Conics: The Most Difficult Math Problem I Ever Solved (CodeParade, 2020).

Consider the following quadrilateral, completed by extending some of its edges:

It can be regarded as a perfect square, drawn in perspective. The "undistorted" square, along with its two points at infinity, are drawn below:

Then, the question is simply to find a projective transform that takes the undistorted square to the distorted square. This can be done with the following matrix:
$$M= \begin{pmatrix}
k_1a_1 & k_2a_2 & k_3a_3 \\
k_1b_1 & k_2b_2 & k_3b_3 \\
k_1    & k_2    & k_3 
\end{pmatrix}$$
where $\vec k = (\vec a - a_0 \vec 1) \times (\vec b - b_0 \vec 1)$, and $\vec a = (a_1, a_2, a_3), \vec b = (b_1, b_2, b_3), \vec 1 = (1, 1, 1)$.

The matrix equation for the inscribed circle is in projective coordinates
$$ \begin{pmatrix} x & y & z \end{pmatrix}
\begin{pmatrix} 1 & 0 & -1/2 \\ 0 & 1 & -1/2 \\ -1/2 & -1/2 & 1/4 \end{pmatrix}
\begin{pmatrix} x \\ y \\ z \end{pmatrix} = 0$$
so after the projective transform, it becomes
$$ \begin{pmatrix} x & y & z \end{pmatrix} N^T
\begin{pmatrix} 1 & 0 & -1/2 \\ 0 & 1 & -1/2 \\ -1/2 & -1/2 & 1/4 \end{pmatrix} N
\begin{pmatrix} x \\ y \\ z \end{pmatrix} = 0$$
where $ N = M^{-1}$

Thus, the coefficient matrix for the transformed circle (now an ellipse) is
$$N^T\begin{pmatrix} 1 & 0 & -1/2 \\ 0 & 1 & -1/2 \\ -1/2 & -1/2 & 1/4 \end{pmatrix} N$$
from which the ellipse can easily be drawn.