# Functions that map a quadrilateral to the unit square?

Given some quadrilateral $$Q \subset \mathbb R^2$$ defined by the vertices $$P_i = (x_i,y_i), i=1,2,3,4$$ (you can assume they are in positive orientation), is there a function $$f: \mathbb R^2 \to \mathbb R^2$$ that is particularly easy to compute which satisfies

$$f(Q) \subseteq U \quad \text{ and } \quad f(\mathbb R^2 \setminus Q) \subseteq \mathbb R^2 \setminus U?$$

Here $$U = [0,1]\times[0,1]$$ denotes the unit square (or $$U=[-1,1]\times[-1,1]$$ if you prefer)

My first attempt was using a function $$f(x,y) := (a + bx + cy +dxy, a' + b'x+c'y +d'xy)$$ (known as the 2d perspecitve transformation bilinear interpolation), but determining the coefficients $$a,b,c,\ldots, d'$$ requires inveriting $$4\times 4$$ matrix to solve two linear systems of equations.

EDIT: The actual 2d perspective transformation as described here does only produce the desired result if $$Q$$ is convex, which is not necessarily the case.

• What you’ve written down is not a planar perspective transformation because of the terms in $xy$. It looks more akin to bilinear interpolation. Anyway, computing a planar perspective transformation between two quadrilaterals doesn’t require a $4\times4$ inversion, if for no other reason than that you’re working with $3\times3$ matrices. See this for an easy-to-compute method. – amd Jan 11 at 0:44
• Sorry, I mixed up the terms, (after looking itup: bilinear interpolation is also used in computer graphics for the transformation of e.g. textures as it preserves axis aligned lines, but it is not the exactly the same as the perspective transformation). Thanks for linking this answer. It does look quite simple indeed, but still requires two inverses of 3×3 matrices. I'll have to do the calculation, but I suspect we still get quite large terms to compute. – flawr Jan 11 at 13:07
• Since you’re mapping the the unit square, the “$B$” matrix has a particularly simple form that represents sums and differences of various columns of $A^{−1}$. However, a homography only has the right properties for a convex quad, anyway, as you point out. Note that if you do have a concavity, you need to take a bit of care with vertex choices for a bilinear transformation-based solution. – amd Jan 11 at 20:01
• Hmm... non-convex quads might be a problem for the bilinear approach, too. The inverse of the bilinear map from the unit square to such a quad is undefined for a large region of $\mathbb R^2\setminus Q$ (I’m assuming that the condition in your question was meant to be $f(\mathbb R^2\setminus Q)\subseteq \mathbb R^2\setminus U$). You might need to do something like break $Q$ up into a pair of triangles. – amd Jan 11 at 20:18
• @amd Yep thanks for the correction, I added it. – flawr Jan 11 at 22:03

The first attempt was bound to fail, as the bilinear transformation suggested only perserves axis aligned lines when mapping a rectangle to some general quadrilateral. The inverse does not have the same form but it is not too difficult to compute, but still not particularly easy. Just for completeness sake here is how this works: Let us first focus on the forward mapping:

Forward mapping: We first define $$g: \mathbb R^2 \to \mathbb R^2$$ with $$g|_U : U \to Q$$ with the bilinear ansatz. It turns out we can easily write this as

\begin{align} (u,v) := g(x,y) &= (1-x)(1-y)P_0 + x(1-y)P_1 + x y P_2 + (1-x)yP_3 \\ &= \underbrace{P_0}_{=(a,a')} + \underbrace{[-P_0+P_1]}_{=(b,b')}x + \underbrace{[-P_0+P_3]}_{=(c,c')}y + \underbrace{[P_0-P_1+P_2-P_3]}_{=(d,d')}xy \end{align}

and we immediately find the coefficients $$a,b,\ldots,d'$$.

Inverse mapping: The inverse $$f:\mathbb R^2 \to \mathbb R^2$$ with $$f|_Q : Q \to U$$ of $$g$$ can then be computed by isolating $$x$$ from the equation for $$u$$:

$$x = \frac{u-a-cy}{b+dy}$$

Pluggin this back into the quation for $$v$$ and multiplying by the denominator we get

$$v(b+dy) = a' (b+dy) + b'(u-a-cy) + c' y(b+dy) + d' (u-a-cy)$$

which we can rearrange as

$$Ay^2+By+C = 0$$

with

\begin{align*} A &= b(a'-v)+b'(u-a) \\ B &= d(a'-v)+d'(u-a) + bc' - cb' \\ C &= dc'-cd' \end{align*}.

This lets us compute $$y$$ and therefore $$x$$ for some given $$(u,v)$$.

EDIT: As @amd commented, this doesn't work for nonconvex quadrilaterals, here the plot of such an example:

• Very nice. This has some problems when $Q$ is not convex, though. No matter how you order the vertices, the forward map sends some points of $U$ to points outside of $Q$. – amd Jan 11 at 20:44
• @amd You're totally right, I made a plot that visualizes this problem! – flawr Jan 11 at 22:02

I expect that we have rather different thresholds for “easy to compute.” A small matrix inversion or two might get a bit tedious, but it’s not very difficult. Both a bilinear map and a planar homography are fairly easy to compute when one of the quadrilaterals is the unit square. Unfortunately, they only meet the criteria of the question when $$Q$$ is convex.

To handle the general case, I think it’s simplest to split the source quad $$Q$$ into a pair of triangles, and correspondingly split $$U$$ along the diagonal from the origin to $$(1,1)$$. It’s then a simple matter to construct a pair of affine transformations of these $$Q$$-halves onto the $$U$$-halves such that they agree along the diagonal.

W.l.o.g. let’s say that we split $$Q$$ along $$P_1P_3$$, so that we have $$L=\triangle{P_1P_2P_3}$$ and $$R=\triangle{P_1P_3P_4}$$. Let $$H$$ be the half-plane, including the boundary line $$P_1P_3$$, in which $$P_2$$ lies and construct $$\mathtt M_L = \begin{bmatrix}0&1&1 \\ 0&0&1 \\ 1&1&1\end{bmatrix}\begin{bmatrix}P_1&P_2&P_3\\1&1&1\end{bmatrix}^{-1} \\ \mathtt M_R = \begin{bmatrix}0&1&0 \\ 0&1&1 \\ 1&1&1\end{bmatrix}\begin{bmatrix}P_1&P_3&P_4\\1&1&1\end{bmatrix}^{-1}.$$ The matrix multiplications here are almost trivial to compute, as the resulting columns are simple sums of the columns of the right-hand matrix. Working in homogeneous coordinates, the map is then defined piecewise: $$f(P) = \begin{cases}\mathtt M_L P &\text{if }P\in H \\[2ex] \mathtt M_R P &\text{otherwise.}\end{cases}$$ Use whatever method is convenient to determine the half-plane in which $$P$$ lies, such as checking the sign of $$P\cdot(P_3\times P_1)$$ (again in homogeneous coordinates).

Or, if you prefer to work with algebraic expressions, construct the inverse map: $$f^{-1}(u,v) = \begin{cases} P_1+u(P_2-P_1)+v(P_3-P_2) &\text{if }u\le v \\[2ex] P_1+u(P_3-P_4)+v(P_4-P_1) &\text{otherwise}\end{cases}$$ and invert those formulas.