Linear interpolation of three 3D points in 3D space

Question

I have three 3D points like p1(x1,y1,z1), p2(x2,y2,z2), p3(x3,y3,z3). I have another point, but I know only x, y value of that point like p4(x4,y4,Z), in which Z is the value I like to compute.

I am sure p4(x4,y4) point is inside a triangle formed by p1(x1,y1), p2(x2,y2), p3(x3,y3) by checking with delaunay triangulation approach. How can I compute Z value of point p4? I like to implement it in C programming. Actually I am trying to implement griddata in MATLAB.

Thanks

U need to provide more information, p1(x1,y1),p2(x2,y2),p3(x3,y3) are they the same points as the 3D points ? Or they are the p1,p2,p3 projections to xy plane ? — tomriddle_1234, Sep 12 '13 at 04:39
z can have multiple solutions. For example, say the known points are (0, 10, 0,) (0, 10, 10) (0, 0, 10) and you know the 4th point as (0, 5, z ), then z can have a set of values from 5 to 10. — Abhishek Bansal, Sep 12 '13 at 05:48
This is looking for the intersection of a plane and a line. As @AbhishekBansal notices, that may have an infinite number of solutions, or it may have zero when the two are parallel. — MSalters, Sep 12 '13 at 07:43

score 7 · Answer 1 · edited Sep 12 '13 at 06:59

7

p1, p2, p3 define a plane. You can represent it by a point and a normal. For instance, P=p1, N=(p2-P) x (p3-P) (that is, N = cross product of p1p2 and p1p3).

Now for p4 to be in the same plane, it satisfies the plane equation:

  (p4-P) · N = 0  %// dot product
⇒ (x4-x1)*N.x + (y4-y1)*N.y + (z4-z1)*N.z = 0

Re-arranging:

z4 = z1 - ((x4-x1)*N.x + (y4-y1)*N.y)/ N.z

No linear system to solve, you just need a cross product.

edited Sep 12 '13 at 06:59

Rody Oldenhuis

37,726
7
50
96

answered Sep 12 '13 at 06:35

DanielKO

4,422
19
29

+1, but write out the solution to the linear system and you'll find the same outcome. – Rody Oldenhuis Sep 12 '13 at 06:57

MBo · Accepted Answer · 2013-09-12T06:15:59.860

5

You can express P4 coordinates in the P1P2P3 vector basis.

x4 = x1 + A * (x2 - x1) + B * (x3 - x1)
y4 = y1 + A * (y2 - y1) + B * (y3 - y1)

This is easy-to-solve linear equation system. You have to find A and B coefficients, then use them to calculate z-coordinate

z4 = z1 + A * (z2 - z1) + B * (z3 - z1)

edited Sep 12 '13 at 06:15

answered Sep 12 '13 at 04:28

MBo

77,366
5
53
86

Yeah it is same as Barycentric Coordinates interpolation. – batuman Sep 12 '13 at 06:31
@MBo: Perhaps due to the gratuitous complexity of solving a 2x2 system when the value can be computed directly? – DanielKO Sep 12 '13 at 06:42
@DanielKO This system may be solved with pen and paper or with Matlab – MBo Sep 12 '13 at 06:59
2

@DanielKO: that's not reason to downvote; it's a valid solution, that I would say has equivalent complexity to your approach. – Rody Oldenhuis Sep 12 '13 at 07:02
For the record, I didn't downvote it, as it's a valid answer; it's just a bit indirect. – DanielKO Sep 12 '13 at 07:56

score 1 · Answer 3 · answered Sep 12 '13 at 06:52

This is to support both MBo's and Konstantin's answers. Please don't accept this question, but one of the others.

This is how you would implement a solution in MATLAB:

%// Your known 3 points
p1 = [ 1 10  0]';
p2 = [-1 10 10]';
p3 = [ 0  0 10]';

%// your 4th target point
p4 = [0 5  NaN]';

%// Difference matrix/vector
A = [p2-p1  p3-p1];
b = p4-p1;

%// Compute solution
p4(end) = p1(end) + A(3,:)*(A(1:2,:)\b(1:2));

Now, in C++, the mere fact of including the relevant eigen libraries blows up the executable size rather spectacularly. What eigen is capable of is complete overkill for this simple 2x2 system.

So I wouldn't go as far as resort to eigen, unless you have tons of other linear algebra things to do. It is a simple 2x2 system, which is easy enough to solve by hand.

Just KISS it; see DanielKO's answer :)

score 0 · Answer 4 · edited Jan 07 '14 at 12:37

0

The mathematical problem here is to solve the following system of equations

p1 + a * (p2 - p1) + b * (p3 - p1) = (x4, y4, z4)

or equivalently

 (x1,y1,z1) + a * (x2-x1, y2-y1, z2-z1) + b * (x3-x1, y3-y1, z3-z1) = (x4, y4, z4)

for a, b, and z4.

To solve it in C/C++, you could either implement the Gauss algo (see also the Numerical Recipes book, it is available online), or use Linear Algebra libraries, such as Eigen, or others.

Remark: the approach is the same regardless if the point (x4, y4) lies within the triangle (x1, y1), (x2, y2), (x3, y3), or not.

edited Jan 07 '14 at 12:37

BenMorel

34,448
50
182
322

answered Sep 12 '13 at 05:43

Konstantin

2,451
1
24
26

thanks. But if (x4,y4) lies inside the three points, the approach is the same. Is it? – batuman Sep 12 '13 at 06:13
@batuman: indeed. If it were to lie *outside* the triangle, it would be called *extrapolation* :) – Rody Oldenhuis Sep 12 '13 at 06:25
Solvable, yes, but that does not mean it always has a solution. – MSalters Sep 12 '13 at 07:45
@MSalters: I don't see the difference between "is solvable" and "has a solution", but you are right in that it does not always have a solution. I edited the answer. – Konstantin Sep 12 '13 at 07:57
@Konstantin: Easiest to see with the roots a second degree polynomial. Always solvable, possibly one solution, but can also have 0 or 2 solutions. – MSalters Sep 12 '13 at 08:15
@MSalters: Still don't get what you mean by "solvable". How can an equation be always solvable if it can have 0 solutions? What is your definition of "solvable"? – Konstantin Sep 12 '13 at 08:22
@Konstantin: An equation is solvable if there exists an algorithm to produce all its solutions, or prove that there are no solutions. `erf(x)=1/3` is not solvable, even though a solution exists. `x*x = -1` is solvable, precisely because we can prove no solution exists. – MSalters Sep 12 '13 at 08:30

Linear interpolation of three 3D points in 3D space

4 Answers4

Linked