Inverse Interpolation for a Quadrilateral

Robert P. Munafo, 2023 Feb 20

This is the inverse of the calculation of two-dimensional interpolation described separately in that article. As in that article the coordinates of the quadrilateral's vertices are known. Here we are given the coordinates of the interpolated point (v, w) and we are trying to determine what parameters p and q were used.

Given a quadrilateral defined by four vertices with coordinates (a, b), (c, d), (e, f), and (g, h), and the coordinates (v, w) of a location within the quadrilateral, we wish to find the parameters p and q that would give (v, w) when using the two-dimensional interpolation calculation. For normal (non-self-crossing) quadrilaterals, if (v, w) is inside the quadrilateral then p and q will be in the range (0..1).

We know that the twelve variables are related by the equations

v = (a(1-p)+c p)(1-q) + (e(1-p)+g p)q (1)
w = (b(1-p)+d p)(1-q) + (f(1-p)+h p)q (2)

and all of the values are known except for p and q. Therefore we can transform the equations (1) and (2) into two equations in p and q:

a - ap + cp - aq + apq - cpq + eq - epq + gpq = v (3)
b - bp + dp - bq + bpq - dpq + fq - fpq + hpq = w (4)

Grouping the p's and q's:

(c-a)p + (e-a)q + (a-c-e+g)pq = v - a (5)
(d-b)p + (f-b)q + (b-d-f+h)pq = w - b (6)

Alter the first equation so that it expresses p in terms of all the rest:

(c-a)p + (a-c-e+g)pq = v-a + (a-e)q     (7)
(c-a + (a-c-e+g)q) p = v-a + (a-e)q     (8)
p = (v-a + (a-e)q)/(c-a + (a-c-e+g)q)     (9)

For simplicity we define new values r, s, t, and u to stand for the combinations of the knowns: use r=v-a, s=a-e, etc. so that (9) becomes

p = (r+sq)/(t+uq) (10)

Similarly for simplicity we define new values j, x, y, and z to stand for the combinations of known values in (6), which becomes

jp + xq + ypq = z (11)

Now use (10) to substitute for p in (11) to get

j(r+sq)/(t+uq) + xq + y(r+sq)q/(t+uq) = z (12)

Multiply everything by (t+uq) to get

j(r+sq) + x(t+uq)q + y(r+sq)q = z(t+uq) (13)

We can regroup to express it as a second-degree polynomial in q, because everything else is a combination of known constants:

jr + jsq + xtq + xuq² + yrq + ysq² = zt + zuq
(xu+ys)q² + (js+xt+yr-zu)q + (jr-zt) = 0

Use the quadratic formula to solve for q, knowing that we want a value between 0 and 1. Once that is found, substitute it back into (10) to get p.

Program

The following program sets up an example and prints the calculations as it goes; you can use the printed values to verify that your implementation is working. The example values are the same as in the two-dimensional interpolation program.

(NOTE: The variable names a, b, ..., h are similar to those in the program in the Jordan Curve Theorem article, but not in the same order; here we have assigned them to the vertices in a "left-to-right, top-to-bottom" order.)

// Perform inverse interpolation for a quadrilateral // // p=?? (9,9) // (5.8,8.2) // (1,7) // (5.96,6.52) q=?? // // // // // (9,1) // (6.6,-0.2) // (3,-2) a = 1 b = 7 c = 9 d = 9 e = 3 f = -2 g = 9 h = 1 v = 5.96 w = 6.52 // compute r,s,t,u so that // p = (v-a + (a-e)q)/(c-a + (a-c-e+g)q) // becomes // p = (r+sq)/(t+uq) r = v-a s = a-e t = c-a u = a-c-e+g print "should be equal:",(v-a+(a-e)*q)/(c-a+(a-c-e+g)*q), (r+s*q)/(t+u*q) // compute j,x,y,z so that // (d-b)p + (f-b)q + (b-d-f+h)pq = w-b // becomes // jp + xq + ypq = z j = d-b x = f-b y = b-d-f+h z = w-b print "should be equal:", (d-b)*p + (f-b)*q + (b-d-f+h)*p*q, w-b print " and equal to:", j*p + x*q + y*p*q, z // substitute "(r+sq)/(t+uq)" for p in "jp+xq+ypq=z" to get an equation // in q: j(r+sq)/(t+uq)+xq+y(r+sq)q/(t+uq) = z print " and equal to:", j*(r+s*q)/(t+u*q)+x*q+y*(r+s*q)*q/(t+u*q) // multiply both sides of the equation by (t+uq) to get // j(r+sq) + x(t+uq)q + y(r+sq)q = z(t+uq) print "should be equal:", j*(r+s*q) + x*(t+u*q)*q + y*(r+sq)*q, z*(t+u*q) // express it as a 2nd degree polynomial in q print "should be 0:", (xu+ys)q^2 + (js+xt+yr-zu)q + (jr-zt) // compute quadratic equation coefficients A,B,C a = x*u+y*s b = j*s+x*t+y*r-z*u c = j*r-z*t print "should be 0:", a*q*q + b*q + c // compute determinant D d = b*b - 4*a*c print "determinant B^2-4AC = ", d if (qd < 0) then print "No solution - determinant is negative" i2i_error = -1 end program end if // try first root q = (-b+sqrt(d))/(2*a) print "root 1: q = ", q if ((q >= 0) and (q <= 1)) then goto GOT_Q // try other root q = (-b-sqrt(d))/(2*a) print "root 2: q = ", q if ((q >= 0) and (q <= 1)) then goto GOT_Q print "neither root is in range" i2i_error = -2 end program GOT_q: // if we get here the is a solution for Q print "got q = ", q if (t+uq) == 0 then print "cannot solve for P (division by zero) i2i_error = -3 end program end if // compute P using: p = (r+sq)/(t+uq) p = (r+sq)/(t+uq) print "got p = ", p // indicate no error i2i_error = 0

revisions: 20230219 initial description; 20230220 complete working program; 20230221 add note about vertex ordering in the sample program; 20230225 fix formatting errors

Mu-ency main page — index — recent changes — DEMZ

This page was written in the "embarrassingly readable" markup language RHTF, and was last updated on 2023 Feb 26.

s.27