### Intro

You have probably done bilinear interpolation in a quadrilateral more than once already. Say you have a random quadrilateral with vertices

**A**,**B**,**C**and**D**, in any number of dimensions. Then any point**X**inside the quadrilateral can be defined by a couple of coordinates*u*and*v*such that the first interpolates across any two opposite edges of the quad to produce two points**P**and**Q**, and such that the second one linearly interpolates again across the line that joins**P**and**Q**. The nice thing is that the same rules that are used to arrive to**X**can be used for the any data stored in the vertices (say, texture coordinates or colors) so that one gets the bilinearly interpolated data values at the point**X**.### The math

The formula that tells us how to get

**X**is very simple:

**P**is the linear interpolation of

**A**and

**B**in

*u*: P = A + (B-A)⋅u

**Q**is the linear interpolation of

**D**and

**C**in

*u*: Q = D + (C-D)⋅u

**X**is the linear interpolation of

**P**and

**Q**in

*v*: X = P + (Q-P)⋅v

therefore

**X**(u,v) =

**A**+ (

**B**-

**A**)⋅

*u*+ (

**D**-

**A**)⋅

*v*+ (

**A**-

**B**+

**C**-

**D**)⋅

*u*⋅

*v*

with

*u*and

*v*running in the [0..1] interval obviously. So far so good. This is old news, really. However, what happens if we now had our quad, we knew the point

**X**and wanted to extract what are the bilinear parameters

*u*and

*v*?

Well, the easy answer is that we simply have to solve this equation and isolate our parameters

*u*and

*v*. "One equation and two unknowns!", you might protest. Remember that the equation above actually applies to all the coordinates of the system, meaning that the vertices

**A**to

**D**are in fact pairs or triplets of numbers (or vectors in general), so in fact we have lots of equations, as many as dimensions we are working on, for only two unknowns. So the systems is very solvable, and the easy answer happens to be the right one.

To make things easier, let's introduce some intermediate variables, such that:

**E**=

**B**-

**A**

**F**=

**D**-

**A**

**G**=

**A**-

**B**+

**C**-

**D**

**H**=

**X**-

**A**

so that the bilinear interpolation equation becomes

**H**=

**E**⋅

*u*+

**F**⋅

*v*+

**G**⋅

*u*⋅

*v*

Now, we choose any two coordinates for our computations. Probably the best choice is to use the 2D plane for which our quad has bigger surface area under orthographic projection. Say we choose axes

*i*and

*j*anyways ( {i,j} can be {x,y}, {z,x}, {y,z}, ...). Then

*u*= (Hi - Fi⋅v)/(Ei + Gi⋅v)

and therefore we can plug it back into the equation to get

Hj⋅Ei + Hj⋅Gi⋅v = Ej⋅Hi - Ej⋅Fi⋅v + Fj⋅Ei⋅v + Fj⋅Gi⋅

*v*

^{2}+ Gj⋅Hi⋅v - Gj⋅Fi⋅

*v*

^{2}

This is a quadratic equation with coefficients:

k2 = Gi*Fj - Gj*Fi

k1 = Ei*Fj - Ej*Fi + Hi*Gj - Hj*Gi

k0 = Hi*Ej - Hj*Ei

which are pretty much some surface areas (kind of 2D cross products) of different edges. The solution is the usual formula

The discriminant k1

^{2}-4⋅k0⋅k2 is always positive if

**X**is inside the quadrilateral

**ABCD**. Then, one has to take the positive or negative square root such that

*v*is in the interval [0..1]. Once we have

*v*,

*u*can be computed with the expression used in the first variable isolation step.

Take must be taken to the cases where the quad is actual a perfect rectangle. In that case

**AB**and

**CD**are parallel, and therefore

**G**= 0. Then k2 = 0 and the equation is basically not a quadratic equation anymore but a simple linear equation of the form k0 + k1⋅

*v*= 0 so that

*v*= -k0/k1

### The code

Here goes some code for the maths above, specialized for 2d:

float cross( in vec2 a, in vec2 b ) { return a.x*b.y - a.y*b.x; }
vec2 invBilinear( in vec2 p, in vec2 a, in vec2 b, in vec2 c, in vec2 d )
{
vec2 e = b-a;
vec2 f = d-a;
vec2 g = a-b+c-d;
vec2 h = p-a;
float k2 = cross( g, f );
float k1 = cross( e, f ) + cross( h, g );
float k0 = cross( h, e );
float w = k1*k1 - 4.0*k0*k2;
if( w<0.0 ) return vec2(-1.0);
w = sqrt( w );
float v1 = (-k1 - w)/(2.0*k2);
float u1 = (h.x - f.x*v1)/(e.x + g.x*v1);
float v2 = (-k1 + w)/(2.0*k2);
float u2 = (h.x - f.x*v2)/(e.x + g.x*v2);
float u = u1;
float v = v1;
if( v<0.0 || v>1.0 || u<0.0 || u>1.0 ) { u=u2; v=v2; }
if( v<0.0 || v>1.0 || u<0.0 || u>1.0 ) { u=-1.0; v=-1.0; }
return vec2( u, v );
}

And here's a realtime running implementation (click on the title to jump to the source code):