These are the pictures I took:
To recover an equation in least squares format, I did the following: $$p' = Hp$$ $$\begin{bmatrix} wx' \\ wy' \\ w \end{bmatrix} = \begin{bmatrix} a & b & c \\ d & e & f \\ g & h & 1 \end{bmatrix} \cdot \begin{bmatrix} x \\ y \\ 1 \end{bmatrix}$$ From this, we can recover two equations: $$\begin{bmatrix} a \\ b \\ c \end{bmatrix} \cdot \begin{bmatrix} x \\ y \\ 1 \end{bmatrix} = \begin{bmatrix} g \\ h \\ 1 \end{bmatrix} \cdot \begin{bmatrix} x \\ y \\ 1 \end{bmatrix} x'$$ $$\begin{bmatrix} d \\ e \\ f \end{bmatrix} \cdot \begin{bmatrix} x \\ y \\ 1 \end{bmatrix} = \begin{bmatrix} g \\ h \\ 1 \end{bmatrix} \cdot \begin{bmatrix} x \\ y \\ 1 \end{bmatrix} y'$$ Formatting this as \(Ax = b\), where \(x = \begin{bmatrix} a & b & c & d & e & f & g & h\end{bmatrix}\), $$\begin{bmatrix} 1 & 1 & 1 & 0 & 0 & 0 & xx' & yx' \\ 0 & 0 & 0 & 1 & 1 & 1 & xy' & yy' \\ \end{bmatrix} \cdot \begin{bmatrix} a \\ b \\ c \\ d \\ e \\ f \\ g \\ h \end{bmatrix} = \begin{bmatrix} x' \\ y' \end{bmatrix}$$ After stacking these for all of the selected point correspondences, I solved the linear equation with np.linalg.lstsq.
Something's not quite correct with my image warping. It seems to work okay for simple transformations such as an identity, shift, or stretch. However, the recovered homographies for certain situations (e.g. trapezoid->square) cause the image to invert in some areas, and other strange behavior. I'm unclear whether this is due to incorrect homography recovering or incorrect warping procedure, but the recovered homographies seem to at least be accurate for the specified correspondence points.
With whatever errors are going on with image warping, the rectification does not look very good :(
Given that warping/rectification aren't fully working yet, I didn't see a point to mosaicing just yet.