Search code examples
algorithminterpolation

Bi-linear interpolation on a trapezoid

I am trying to interpolate data based on an input forming a trapezoid. I found this solution by Nico Schertler, which outlines how to interpolate for any four non-aligned input points. However, I am unable to get it to generalize to trapezoids.

The interpolated value is found through:

interpolated value = (1 - alpha) * ((1 - beta) * p1 + beta * p2) + alpha * ((1 - beta) * p3 + beta * p4)

Where

alpha = -(b e - a f + d g - c h + sqrt(-4 (c e - a g) (d f - b h) +
        (b e - a f + d g - c h)^2))/(2 c e - 2 a g)    
beta  = (b e - a f - d g + c h + sqrt(-4 (c e - a g) (d f - b h) + 
        (b e - a f + d g - c h)^2))/(2 c f - 2 b g)

or

alpha = (-b e + a f - d g + c h + sqrt(-4 (c e - a g) (d f - b h) + 
        (b e - a f + d g - c h)^2))/(2 c e - 2 a g)
beta  = -((-b e + a f + d g - c h + sqrt(-4 (c e - a g) (d f - b h) + 
        (b e - a f + d g - c h)^2))/( 2 c f - 2 b g))

The coefficients are given by:

a = -p1.x + p3.x
b = -p1.x + p2.x
c = p1.x - p2.x - p3.x + p4.x
d = interpolated_point.x - p1.x
e = -p1.y + p3.y
f = -p1.y + p2.y
g = p1.y - p2.y - p3.y + p4.y
h = interpolated_point.y - p1.y

Where p1-p4 are coordinates of input vertices of the polygon. The issue occurs when p1.x=p3.x and p2.x=p4.x but p1.y!=p2.y!=p3.y!=p4.y - which is the case for a trapezoid. In this case, a=c=0 and the alpha coefficient becomes becomes infinite upon division by 0.

So my question is: is there any way to modify (simplify?) Nico's method for it to be applicable to a trapezoid.

Solution

  • There is only one unique solution for a trapezoid, which implies that the equations for the parametric coordinates are linear. Therefore the general solution for an arbitrary quadrilateral is undefined since the equations it was derived from are quadratic.

    This can of course be proven rigorously, but can also be easily illustrated with two diagrams:

    enter image description here

    enter image description here

    In both configurations (which are equivalent if the parameters are swapped), the blue lines representing different values of alpha never intersect inside the trapezoid itself, which means there cannot be multiple solutions.

    The special case solution for trapezoids can be derived from the original expressions for the interpolated coordinate R:

    R = (1 - β)*[(1 - α)*P1 + α*P3] + β*[(1 - α)*P2 + α*P4]

  = (1 - β)*[P1 + α*(P3 - P1)] + β*[P2 + α*(P4 - P2)]

  = P1 + β*(P2 - P1) + α*(P3 - P1) + αβ*(P4 - P3 - P2 + P1)

  = P1 + β*E12 + α*E13 + αβ*(E24 - E13)

where the "edge vector" Exy = Py - Px

    Moving P1 to the LHS, and taking the cross-product of both sides with edges 1-2 and 1-3:

    (R - P1)xE12 = β*E12xE12 + α*E13xE12 + αβ*(E34xE12 - E12xE12)
------------     -------     -------       -------   -------
      A            = 0          C             E        = 0
   
(R - P1)xE13 = β*E12xE13 + α*E13xE13 + αβ*(E24xE13 - E13xE13)
------------     -------     -------       -------   -------
      B             D          = 0            F        = 0

    Recalling that the cross-product of parallel vectors equals 0. Arriving at a pair of simultaneous quadratic equations:

    A = α*C + αβ*E
B = β*D + αβ*F

note that D = -C

    In the first diagram, E13 and E24 are parallel, which means F = 0:

    β = B/D
α = A / (C + E*B/D) = A*C / (C*C - E*B)

    ... and ditto for the second.