Relative position of a point within a quadrilateral

Question

I am trying to find the easiest way to determine a relative position of a point within a quadrilateral. The known are (see figure) the positions of points 1, 2, 3, 4 and 5 in the xy-coordinate system: x1, y1, x2, y2, x3, y3, x4, y4, x5, y5.

Also known are the positions of points 1, 2, 3, and 4 in the ξ-η coordinate systems (see figure).

From this data, I want to determine what are the ξ and η for point 5.

Results

Thank you to all who anwsered! I find the solution by @dbc and @agentp similar. Also I find this solution better than the perspective transformation solution by @MBo, since I do not have to compute the inverse of a matrix (Ax=B --> x=inv(A)*B).

I get the following result for:

u = 0.5 * (ξ + 1)
v = 0.5 * (η + 1)

In my case all points are within the rectangle, therefore u>0 and v>0.

Answer 1

What you have here is a 2d bilinear blended surface. For simplicity, let's change its coordinates to range from zero to one:

u = 0.5 * (ξ + 1)
v = 0.5 * (η + 1)

In that case, the surface evaluator can be expressed as

F(u, v) = P1 + u * (P2 - P1) + v * ((P4 + u * (P3 - P4)) - (P1 + u * (P2 - P1)))

I.e., for a given u, construct a line passing through the following two points:

Pv0 = P1 + u * (P2 - P1);
Pv1 = P4 + u * (P3 - P4);

then interpolate between then for given v

F(u, v) = Pv0 + v * (Pv1 - Pv0)

What you seek are values (u,v) such that F(u, v) = P5. This will occur for given u when the line from Pv0 to Pv1 passes through P5, which will occur when P5 - Pv0 is parallel to Pv1 - Pv0 -- i.e. when their 2d cross is zero:

cross2d(P5 - Pv0, Pv1 - Pv0) = 0

⇒

cross2d(P5                 - (P1 + u * (P2 - P1)),  
        P4 + u * (P3 - P4) - (P1 + u * (P2 - P1))) = 0

Now, the 2d cross of two 2d vectors A ⨯ B is given by Ax*By - Ay*Bx, so that equation becomes

(x5 - (x1 + u * (x2 - x1))) * (y4 + u * (y3 - y4) - (y1 + u * (y2 - y1))) - (y5 - (y1 + u * (y2 - y1))) * (x4 + u * (x3 - x4) - (x1 + u * (x2 - x1))) = 0

Expanding this expression out and collecting collecting together terms in u, we get

    u^2 * (x1*y3 - x1*y4 - x2*y3 + x2*y4 + (-x3)*y1 + x3*y2 + x4*y1 - x4*y2)
  + u *   (-x1*y3 + 2*x1*y4 - x1*y5 - x2*y4 + x2*y5 + x3*y1 - x3*y5 - 2*x4*y1 + x4*y2 + x4*y5 + x5*y1 - x5*y2 + x5*y3 - x5*y4)
  +       (-x1*y4 + x1*y5 + x4*y1 - x4*y5 - x5*y1 + x5*y4)
= 0

This is now a quadratic equation over u, and can be solved as such. Note that in cases where the top and bottom edges of your quadrilateral are parallel then the quadratic devolves into a linear equation; your quadratic equation solver must needs handle this.

        double a = (x1 * y3 - x1 * y4 - x2 * y3 + x2 * y4 + (-x3) * y1 + x3 * y2 + x4 * y1 - x4 * y2);
        double b = (-x1 * y3 + 2 * x1 * y4 - x1 * y5 - x2 * y4 + x2 * y5 + x3 * y1 - x3 * y5 - 2 * x4 * y1 + x4 * y2 + x4 * y5 + x5 * y1 - x5 * y2 + x5 * y3 - x5 * y4);
        double c = (-x1 * y4 + x1 * y5 + x4 * y1 - x4 * y5 - x5 * y1 + x5 * y4);

        double[] solutions = Quadratic.Solve(a, b, c);

There may be more than one solution. There might also be no solutions for a degenerate quadrilateral.

Having solved for value(s) of u, finding the equivalent v is straightforward. Given points

Pv0 = P1 + u * (P2 - P1);
Pv1 = P4 + u * (P3 - P4);

you seek v such that

v * (Pv1 - Pv0) = P5 - Pv0;

Pick the coordinate index 0 or 1 such that |(Pv1 - Pv0)[index]| is maximized. (If both coordinates are almost zero, then give up -- there's no solution for this specific u. Then set

v = (P5 - Pv0)[index] / (Pv1 - Pv0)[index];

Finally, if you have more that one solution, prefer a solution inside the [u, v] boundaries of the blend. Then finally set

ξ = 2 * u - 1;
η = 2 * v - 1;

Answer 2

This looks like a standard finite element parameterization (The question doesn't specify a particular mapping, but I imagine someone might be interested in this specific case)

 {x, y}  == (
             (1 - eta) (1 - ci) {p1x, p1y}  +
             (1 - eta) (1 + ci) {p2x, p2y}   + 
             (1 + eta) (1 + ci) {p3x, p3y} +
             (1 + eta) (1 - ci) {p4x, p4y} )/4

This can be solved in closed form for {eta,ci}, but the expression is pretty unwieldy to post.

In practice, compute these constants:

 ax = p1x + p2x + p3x + p4x
 bx = p1x - p2x - p3x + p4x
 cx = p1x + p2x - p3x - p4x
 dx = p1x - p2x + p3x - p4x
 ay = p1y + p2y + p3y + p4y
 by = p1y - p2y - p3y + p4y
 cy = p1y + p2y - p3y - p4y;
 dy = p1y - p2y + p3y - p4y;

Solve this quadratic for eta :

 (ax by - bx ay) - 4 (by  x - bx y) +
  eta (dx ay - cx by + bx  cy - ax dy + 4 (x dy - dx y)) +
  eta^2 (cx dy - dx cy) == 0

then get ci as:

 ci = ((-ax + eta cx + 4 x)/(-bx + eta dx))

If the polygon is not too distorted just one of the solutions will satisfy -1<eta<1 and -1<ci<1

Answer 3

Referring to the self-answer of @blaz (please vote up the answers of @blaze, @dbc and @agentp)

For everybody who is not willing to copy the formulas by hand, here is the formula as C# code:

    double v_sqrt = Math.Sqrt(
        4 * (
        (x3 - x4) * (y1 - y2) - (x1 - x2) * (y3 - y4)) * (x4 * (-1 * y + y1) + x1 * (y - y4) + x * (-1 * y1 + y4)) +
        Math.Pow(
        (x3 * y - x4 * y - x3 * y1 + 2 * x4 * y1 - x4 * y2 + x1 * (y + y3 - 2 * y4) + x2 * (-1 * y + y4) + x * (-1 * y1 + y2 - y3 + y4))
        , 2)
        );

    double u_sqrt = Math.Sqrt(
            4 * ((x3 - x4) * (y1 - y2) - (x1 - x2) * (y3 - y4))
            * (
                x4 * (-1 * y + y1) + x1 * (y - y4) + x * (-1 * y1 + y4)
            ) +
            Math.Pow(
            (x3 * y - x4 * y - x3 * y1 + 2 * x4 * y1 - x4 * y2 + x1 * (y + y3 - 2 * y4) + x2 * (-1 * y + y4) + x * (-1 * y1 + y2 - y3 + y4))
            , 2)
        );

    double k = 1 / (2 * ((x3 - x4) * (y1 - y2) - (x1 - x2) * (y3 - y4)));
    double l = 1 / (2 * ((x1 - x4) * (y2 - y3) - (x2 - x3) * (y1 - y4)));

    ///////////////////////////////////////////////////////////////////////////////////////////////
    double v1 = l *
        (x2 * y - x3 * y + x4 * y + x * y1 - 2 * x2 * y1 + x3 * y1 - x * y2 - x4 * y2 + x * y3 - x1 * (y - 2 * y2 + y3) - x * y4 + x2 * y4 +
        v_sqrt);

    ///////////////////////////////////////////////////////////////////////////////////////////////
    double u1 = -1 * k *
        (-x2 * y + x3 * y - x * y1 - x3 * y1 + 2 * x4 * y1 + x * y2 - x4 * y2 - x * y3 + x1 * (y + y3 - 2 * y4) + x * y4 + x2 * y4 +
        u_sqrt);

    double v2 = -1 * l *
        (x1 * y + x3 * y - x4 * y - x * y1 - 2 * x3 * y1 + x * y2 - -2 * x1 * y2 + x4 * y2 - x * y3 + x1 * y3 + x * y4 - x2 * (y - 2 * y1 + y4) +
        v_sqrt);

    /////////////////////////////////////////////////////////////////////////////////////////////////
    double u2 = k *
        (x2 * y - x3 * y + x4 * y + x * y1 + x3 * y1 - 2 * x4 * y1 - x * y2 + x4 * y2 + x * y3 - x1 * (y + y3 - 2 * y4) - x * y4 - x2 * y4 +
        u_sqrt);

In most cases it is u1 and v1 so there should not be the need for computing the other ones.

I used it to calibrate the coordinates of a Pegasus Air-Pen device (ultrasonic stylus) on a sheet of paper. It does work best if your coordinates for point 1 to 5 are also >= 0.

Sry for posting this as an answer but it is too long for a comment and I think it is a valuable help for this post as it would be for me.

Answer 4

You need to calculate a matrix of perspective transformation, that maps 4 points of source quadrilateral to 4 points of destination quadrilateral (example) (more mathemathics), then apply this transformation to coordinates of 5th point (multiply matrix by coordinate vector)