# Mapping Irregular Quadrilateral to a Rectangle

I have a camera looking at a computer monitor from varying angles. Since the camera is a grid of pixels, I can define the bounds of the monitor in the camera image as:

I hope that makes sense. What I want to do is come up with an algorithm to translate points within this shape to this:

I have points within the same domain as ABCD, as determined from the camera, but I need to draw these points in the domain of the monitor's resolution.

Does that makes sense? Any ideas?

-
You're asking for an affine transformation? – J. M. Dec 7 '10 at 18:43
Yes, that's quite close to what I'm looking for. I'm looking to map a 4 sided shape with corners ABCD to A'B'C'D'. Is there a way to do this in general, or do i need to figure out the individual transform, rotate, shear, etc. differences between the two shapes? – bufferz Dec 7 '10 at 18:56
I have a hard time understanding what you are asking. It seems you have an arbitrary quadrilateral and you want to map it to a rectangle? The most straight-forward way would be to average the opposite sides and take these to be your widths and heights. Is there something you want to be preserved by this transformation, e.g. area? – WWright Dec 7 '10 at 19:17
Not all rectangles are achievable by pointing a camera at a rectangle. So ignore one of the corner-points and then do a linear transform (along with a translation to match up the origin) to match up the remaining triangle. The fourth point will magically end up where it should be. I would make this an answer, but I can't be bothered to actually work out which matrix I need. – Oscar Cunningham Dec 7 '10 at 19:25
A 2D->2D affine transformation will not work, as they map parallelograms to parallelograms. You can do 3D->2D taking advantage of the depth to account for the non-parallelism (the two sides will be parallel in 3 space). If you know the depth, you can just follow the Wikipedia article and solve for the matrix elements of the transform. – Ross Millikan Dec 7 '10 at 19:25

In general there is no affine transformation that maps an arbitrary quadrangle onto a rectangle. But there is (exactly one) projective transformation $T$ that maps a given quadrangle $(A, B, C, D)$ in the projective plane onto a given quadrangle $(A', B', C' D')$ in the same or another projective plane. This $T$ is ${\it collinear}$, i.e., it maps lines to lines. To do the calculations you have to introduce homogeneous coordinates $(x,y,z)$ such that $D=(0,0,1)$, $C=(1,0,1)$, $A=(0,1,1)$, $B=(1,1,1)$ and similarly for $A'$, $B'$, $C'$, $D'$. With respect to these coordinates the map $T$ is linear and its matrix is the identity matrix.

-
I'm a bit rusty on my matrix math. I'm not looking for a handout algorithm but I'm wondering if there's a wikipedia or similar article regarding your suggestion. I'm trying to end up with a function relating a ABCD and A'B'C'D' directly, so I can implement that into my software project. – bufferz Dec 8 '10 at 2:43
It cannot be done without some matrix manipulation. The following page seems to contain very explicit formulas pertaining to your problem: geometrictools.com/Documentation/PerspectiveMapping.pdf – Christian Blatter Dec 8 '10 at 20:04
The link has changed to geometrictools.com/Documentation/PerspectiveMappings.pdf – VoteCoffee Oct 22 '14 at 16:45

The best solution I've found so far on a forum lost in the sea of forums is to decompose your problem like this :

Here, U and V represent coordinates within the quadrilateral (scaled between 0 and 1).

From $P0$, $P1$, $P2$ & $P3$ we can easily compute the normalized normal vectors $N0$, $N1$, $N2$ & $N3$. Then, it's easy to see that : $$u = \frac{dU0}{dU0 + dU1} = \frac{(P-P0) \cdot N0}{(P-P0).N0 + (P-P2) \cdot N2} \\ v = \frac{dV0}{dV0 + dV1} = \frac{(P-P0) \cdot N1}{(P-P0).N1 + (P-P3) \cdot N3}.$$

This parametrization works like a charm and is really easy to compute within a shader for example. What's tricky is the reverse: finding $P(x,y)$ from $(u,v)$ so here is the result:

$$x = \frac{vKH \cdot uFC - vLI \cdot uEB}{vJG \cdot uEB - vKH \cdot uDA}, \\ y = \frac{vLI \cdot uDA - uFC \cdot vJG}{vJG \cdot uEB - vKH \cdot uDA},$$

where: $$uDA = u \cdot (D-A), \quad uEB = u \cdot (E-B), \quad uFC = u \cdot (F-C), \\ vJG = v \cdot (J-G), \quad vKH = v \cdot (K-H), \quad vJG = v \cdot (J-G),$$

and finally: $$A = N0_x, \qquad \qquad B = N0_y, \quad C = -P0 \cdot N0, \qquad \\ D = N0_x + N2_x, \quad E = N0_y + N2_y, \quad F = -P0 \cdot N0 - P2 \cdot N2, \\ G = N1_x, \qquad \qquad H = N1_y, \quad I = -P0 \cdot N1, \qquad \\ J = N1_x + N3_x, \quad K = N1_y + N3_y, \quad L = -P0 \cdot N1 - P2 \cdot N3.$$

I've been using this successfully for shadow mapping of a deformed camera frustum mapped into a regular square texture and I can assure you it's working great! :D

-
What are $u,v$? – nbubis Mar 21 '14 at 0:21

Try this solution, it worked for me.

-
As far as I can tell the solution you point to fails if you have a parallelogram as an input. For example, try giving it a . It's quite possible I'm not understanding something about the solution but I'd love to know what it is if you can get it to work for parallelograms. – Mattia Nov 29 '11 at 21:27
@Mattia The linked solution describes a method (bilinear interpolation) to map the unit square onto an arbitrary quadrilateral. There are no restrictions on the quadrilateral. This is probably exactly what the OP is looking for. If the inverse of this mapping is required, this is also discussed in the link. – yasmar Feb 1 '12 at 14:09

Here is a solution implemented in VBA, a General Algebraic solution, more general than the augmented 2D affine transformation formulation on Wikipedia.

Function Quad_to_Logical_Cell(Qx() As Double, Qy() As Double, x As Double, y As Double) As Variant

'WJW 7-13-15
'This function performs a coordinate transform from X,Y space to the normalized L,M.
'
'If a point {is within {0,1} on both axes, it is within the transformed unit square.
'Qx,Qy vectors contain the 4 coordinates of the corners - x and y values, respectively, ordered as indicated below:
'
'The unit cell L(l,m) corresponding to Q(x,y) is oriented as:
'L0(x=0,y=0),L1(0,1), L2(1,1), L3(1,0).  The order matters.
'The following represent an algebraic solution to the system:
'l=a1 + b1x + c1y + d1xy
'm=a2 + b2x + c2y + d2xy

Dim L_Out() As Double
ReDim L_Out(2)

ax = (x - Qx(0)) + (Qx(1) - Qx(0)) * (y - Qy(0)) / (Qy(0) - Qy(1))
a3x = (Qx(3) - Qx(0)) + (Qx(1) - Qx(0)) * (Qy(3) - Qy(0)) / (Qy(0) - Qy(1))
a2x = (Qx(2) - Qx(0)) + (Qx(1) - Qx(0)) * (Qy(2) - Qy(0)) / (Qy(0) - Qy(1))
ay = (y - Qy(0)) + (Qy(3) - Qy(0)) * (x - Qx(0)) / (Qx(0) - Qx(3))
a1y = (Qy(1) - Qy(0)) + (Qy(3) - Qy(0)) * (Qx(1) - Qx(0)) / (Qx(0) - Qx(3))
a2y = (Qy(2) - Qy(0)) + (Qy(3) - Qy(0)) * (Qx(2) - Qx(0)) / (Qx(0) - Qx(3))
bx = x * y - Qx(0) * Qy(0) + (Qx(1) * Qy(1) - Qx(0) * Qy(0)) * (y - Qy(0)) / (Qy(0) - Qy(1))
b3x = Qx(3) * Qy(3) - Qx(0) * Qy(0) + (Qx(1) * Qy(1) - Qx(0) * Qy(0)) * (Qy(3) - Qy(0)) / (Qy(0) - Qy(1))
b2x = Qx(2) * Qy(2) - Qx(0) * Qy(0) + (Qx(1) * Qy(1) - Qx(0) * Qy(0)) * (Qy(2) - Qy(0)) / (Qy(0) - Qy(1))
by = x * y - Qx(0) * Qy(0) + (Qx(3) * Qy(3) - Qx(0) * Qy(0)) * (x - Qx(0)) / (Qx(0) - Qx(3))
b1y = Qx(1) * Qy(1) - Qx(0) * Qy(0) + (Qx(3) * Qy(3) - Qx(0) * Qy(0)) * (Qx(1) - Qx(0)) / (Qx(0) - Qx(3))
b2y = Qx(2) * Qy(2) - Qx(0) * Qy(0) + (Qx(3) * Qy(3) - Qx(0) * Qy(0)) * (Qx(2) - Qx(0)) / (Qx(0) - Qx(3))

'Dependent on the way your data is formatted, you may have to swap x and y to get the order right.
'L=L(0) is the x coordinate here (row)
'M=L(1) is the y coordinate here (colum)
L_Out(0) = (ax / a3x) + (1 - a2x / a3x) * (bx - b3x * ax / a3x) / (b2x - b3x * a2x / a3x)
L_Out(1) = (ay / a1y) + (1 - a2y / a1y) * (by - b1y * ay / a1y) / (b2y - b1y * a2y / a1y)

End Function

-

I've been wrestling with a very similar problem in order to determine gradients in an irregular quad grid and needing to map points within arbitrary quadrilaterals to a unit square. In addition, I require the inverse mapping of the x and y axis at the mapped normalized coordinate location back into the quad so I can determine the orientation of the quad grid at that point. i.e. if [x',y'] are the transformed coordinates, I need to be able to do an inverse transform on [0,y'],[1,y'] and [x',0],[x',1]. Here is what I've come up with:

You can divide the quad into two tris, and use affine maps on these individually. This is not difficult. This will create a noticable effect at the division between the two tris, however.

If you want a smooth mapping from a quad to a square (or rectangle), you need to use a non-affine transform such as a projective transform. There are other transforms other than projective that will also work, and also be colinear (preserve straight lines).

If [x1,y1],[x2,y2],[x3,y3],[x4,y4] are the four points in the quad, then the 4x4 matrix B in the the following will yield a mapping into the square (on the RHS) that seems to work and may be easier to compute than the proper 3x3 projective matrix.

%     [x1 y1 x1*y1 1]              [0 0 0 1]
%     [x2 y2 x2*y2 1]   X  B  =    [1 0 0 1]
%     [x3 y3 x3*y3 1]              [0 1 0 1]
%     [x4 y4 x4*y4 1]              [1 1 1 1]


The question I have is that if one does this, and then wants to use the inverse of B to do the inverse transform, how do you calculate the third elements of the location vectors for the orthogonal coordinates. (They are no longer x*y.)

NOTE: If you want to map into any other (arbitrary) quadrilateral (such as a rectangle), then just replace the RHS of what I have above with the new coordinates.

%     [x1 y1 x1*y1 1]              [x1' y1' x1'*y1' 1]
%     [x2 y2 x2*y2 1]   X  B  =    [x2' y2' x2'*y2' 1]
%     [x3 y3 x3*y3 1]              [x3' y3' x3'*y3' 1]
%     [x4 y4 x4*y4 1]              [x4' y4' x4'*y4' 1]

-

Take a look on Gernot Hoffmann's tutorial on image rectification. There are also the special cases (rectangle to quadrilateral) explained.

Another page that helped me discusses 2D perspective transformation (i.e. planar homography).

For a deep understanding of the topic and more numerically stable algorithms, I can only recommend Hartley & Zisserman: Multi-View Geometry in Computer Vision.

-

You may find this example Perl code useful, using the Imager library.

-

You could approach this by using an isoparametric mapping. Say the quadrilateral object is said to be in a $x_{1}-y_1$ coordinate frame, while the rectangle is in a new $x_{2}-y_{2}$ frame. What you can do is find $x_{1}=x_{1}(x_{2},y_{2})$ and $y_{1}=y_{1}(x_{2},y_{2})$ using an interpolation-based mapping.

Say we define each vertex as a 2D vector $\vec{P}_{i}$, we can end up with the following mapping to find an a given $\vec{P}$ as a function of $x_{2}$ and $y_{2}$:

$$\vec{P}(x_{2},y_{2}) = \sum_{i=1}^{4}\vec{P}_{i}h_{i}(x_{2},y_{2})$$

Now, assume point A, $\vec{P}_{1}$, corresponds with $(0,0)$ location, point B, $\vec{P}_{2}$, corresponds with $(width,0)=(w,0)$, etc. With that, we can arrive at the following expressions for $h_{i}$:

$$h_{1}(x_{2},y_2) = \frac{(x_{2}-w)(y_{2}-h)}{wh}$$

$$h_{2}(x_{2},y_2) = \frac{x_{2}(h-y_{2})}{wh}$$

$$h_{3}(x_{2},y_2) = \frac{x_{2}y_{2}}{wh}$$

$$h_{4}(x_{2},y_2) = \frac{(w-x_{2})y_{2}}{wh}$$

Using all this information, you could loop through the rectangle to find the $\vec{P}$ coordinate in the original image that each $(x_2,y_2)$ pixel is associated with, then get the pixel information and drop it into the $(x_2,y_2)$ pixel. As a note, the $h_i$ expressions were found via Lagrangian interpolation procedures.

-