# Bilateral interpolation on irregular grid of points?

#### AlbertHall

Joined Jun 4, 2014
9,869
Can I get useful results by doing bilateral interpolation on an irregular set of points?
Suppose the four points around the point I want to know about form an irregular quadrilateral - no right angles, no parallel sides like the image below. #### MrAl

Joined Jun 17, 2014
7,591
Can I get useful results by doing bilateral interpolation on an irregular set of points?
Suppose the four points around the point I want to know about form an irregular quadrilateral - no right angles, no parallel sides like the image below. Hi,

What is the application, and do you have any pass/fail examples?

#### AlbertHall

Joined Jun 4, 2014
9,869
I have a router which I use to make PCBs. To keep the tool depth correct first I probe the surface of the board to measure the height on a square grid and use bilateral interpolation to find the correct depth for any point on the board. I have written a program which adds this auto-levelling code to the gcode for the board.

My latest idea is to be able to have plated through holes. This entails drilling the holes and any milled slots which need plating then removing the board and tin plating it and replacing the board on the machine. The next operation is to etch the tracks but first the probing must be repeated as the board is now thicker and the plating is not necessarily of even thickness. The problem is that now the points probed have to avoid the drilled holes and milled slots so, generally, a simple square grid will not be possible. I am now trying to come up with a way of modifying my code to handle this situation.

#### WBahn

Joined Mar 31, 2012
25,760
So, if I understand you correctly, you have a grid point somewhere within that quadrilateral for which you want to know the nominal height but you can't actually probe it. So you find four nearby points that you can probe the height of and now want to know if you can use those measurements to determine what the proper height should be at the point on your grid that you can't probe.

If you use three of the points, you can define a plane and then find the height of the point in that plane that is vertically above/below the desired grid point. With four points you can do that with four different planes. You could then average the four to get your best estimate and also make a judgement on how good you think that estimate probably is. If all four lie within some maximum span, then you probably have a good estimate. If one is a significant outlier, then you might throw that one out and possibly note that the something is going on in the vicinity of the measured point that produced the outlier.

#### MrAl

Joined Jun 17, 2014
7,591
I have a router which I use to make PCBs. To keep the tool depth correct first I probe the surface of the board to measure the height on a square grid and use bilateral interpolation to find the correct depth for any point on the board. I have written a program which adds this auto-levelling code to the gcode for the board.

My latest idea is to be able to have plated through holes. This entails drilling the holes and any milled slots which need plating then removing the board and tin plating it and replacing the board on the machine. The next operation is to etch the tracks but first the probing must be repeated as the board is now thicker and the plating is not necessarily of even thickness. The problem is that now the points probed have to avoid the drilled holes and milled slots so, generally, a simple square grid will not be possible. I am now trying to come up with a way of modifying my code to handle this situation.
Hi,

Ok good. I have a couple more questions.

First, is the PCB flat or so nearly flat that you can consider it a single plane?

Second, as WBahn is thinking, do these four points correspond to depths where they may be all different and you are going to use them to determine what depth the actual point is, and the point is in the middle somewhere?

Is the point in the middle somewhere not exactly known or do you know the projected x,y coordinates?

If the board is considered perfectly flat and the projected x,y coordinates are known then i think all you need is three points to determine the plane, and once you do that you should be able to calculate the depth of any point on the entire PC board or at least in the vicinity of the three points. Of course we might also take into account the accuracy of the measurements.

#### AlbertHall

Joined Jun 4, 2014
9,869
If you use three of the points, you can define a plane and then find the height of the point in that plane that is vertically above/below the desired grid point. With four points you can do that with four different planes.
Right. So now I have a new phrase in my vocabuary - barycentric coordinates. The point will not necessarily be in all four of those planes but if any of three weights calculated by this barycentric method are negative that tells you that the point is outside that triangle. That's very convenient.

#### AlbertHall

Joined Jun 4, 2014
9,869
First, is the PCB flat or so nearly flat that you can consider it a single plane?
Not the whole board, no. With the square grid I was using a 1cm grid and assuming the board was flat within each 1cm square. PCBs tend to be curved and it is difficult to get rid of that unless you have a very flat vacuum table that can suck the board flat. I don't have one of those.

Second, as WBahn is thinking, do these four points correspond to depths where they may be all different and you are going to use them to determine what depth the actual point is, and the point is in the middle somewhere?
Yup, exactly.

Is the point in the middle somewhere not exactly known or do you know the projected x,y coordinates?
The point in the middle is the point to be milled. Usually two points within the grid with heights determined, then mill from depth1 at point 1 to depth 2 at point 2. So the cut goes all the way through the copper and as little as possible into the fibreglass as that wrecks tools. Remember the copper is just 1.37mils thick - 0.00137 inches - so depth control needs to be pretty precise.

#### MrAl

Joined Jun 17, 2014
7,591
Hello again,

Ok then once you determine what method you want to use to determine the x,y point you can use this solution:
z=(((b1-1)*c3-b3*c1+b3)*y+((-b2-2)*c3+b3*c2+2*b3)*x+
((b1-1)*c2+(-b2-2)*c1+b2+2*b1)*u3+((1-b1)*c3+b3*c1-b3)*u2+
((b2+2)*c3-b3*c2-2*b3)*u1)/((b1-1)*c2+(-b2-2)*c1+b2+2*b1)

That's a solution for the 'z' coordinate of any point on a plane where the plane is given by three points, a, b, and c, and given the projected x,y coordinate of the point in question.
Also, 'u' is either of the three points you can choose as you see fit.
The depth will of course be relative to the level surface you decide to use for all the measurements.

Just to note:
All of the points a,b,c,u have three dimensions [x,y,z], so for example 'a' might be:
[1,2,3]
and 'b' might be:
[4,2,7]
etc.

Your 1cm constraint could be your max distance between any two points then on the level surface for example.

#### AlbertHall

Joined Jun 4, 2014
9,869
Have I got this right?
b1 is 'x' coordinate of point b.
c3 is 'z' coordinate of point c.

I don't get 'u'. Is it a simple duplicate of any one of the other three points?

#### MrAl

Joined Jun 17, 2014
7,591
Have I got this right?
b1 is 'x' coordinate of point b.
c3 is 'z' coordinate of point c.

I don't get 'u'. Is it a simple duplicate of any one of the other three points?
Hi,

Yes:
a1 is x1, a2 is y1, a3 is z1,
b1 is x2, b2 is y2, b3 is z2,
c1 is x3, c2 is y3, c3 is z3,
u can be a, b, or c, so for example if 'a' then:
u1 is x1, u2 is y1, u3 is z1.

The reason for 'u' is so you can choose either of the three points, and perhaps a fourth point if it is on the plane. Maybe choose 'u' such that it is the most accurately measured, or something like that.

Using that we get:
z=(u1*((y2-y1)*z3+y3*(z1-y3)-y2*z1+y1*y3)
+x*((y1-y2)*z3+y2*z1+y3*(y3-z1)-y1*y3)
+y*((x2-x1)*z3+x3*(z1-y3)-x2*z1+x1*y3)
+u2*((x1-x2)*z3+x2*z1+x3*(y3-z1)-x1*y3)
+u3*((x2-x1)*y3+x1*y2+x3*(y1-y2)-x2*y1))
/((x2-x1)*y3+x1*y2+x3*(y1-y2)-x2*y1)

and choose 'u' as before.

Last edited:

#### AlbertHall

Joined Jun 4, 2014
9,869
The formula in post #8 uses 'x' and 'y' for the location of the point whose height is to be determined, but the post #10 formula doesn't seem to include any reference to the location of that point?

[EDIT] Ignore that. Going blind and mad looking at these nice simple formulas.

#### AlbertHall

Joined Jun 4, 2014
9,869
I have put the post #10 formula and also the barycentric coordinate formula into a spreadsheet, attached below.
I have selected a point near one corner of the triangle so I know the answer should be near the height of that corner.

The barycentric method gives a believable answer, as does MrAl formula unless I set u to match point 2 when it is a long way off.
Presumably I have entered something wrong but I have done the conversion twice and compared them and got the same.

Is there any reason choosing that point should mess up the answer?

#### Attachments

• 37.4 KB Views: 2
• 7.5 KB Views: 1

#### MrAl

Joined Jun 17, 2014
7,591
I have put the post #10 formula and also the barycentric coordinate formula into a spreadsheet, attached below.
I have selected a point near one corner of the triangle so I know the answer should be near the height of that corner.

The barycentric method gives a believable answer, as does MrAl formula unless I set u to match point 2 when it is a long way off.
Presumably I have entered something wrong but I have done the conversion twice and compared them and got the same.

Is there any reason choosing that point should mess up the answer?
Hi,

I cant read your files, can you give me the three points any any other formula you are using in text?
I'll also go over the formula again to make sure it is correct.

[LATER]
Could i have posted the wrong formula?
z=
(((b1-a1)*c3+(a3-b3)*c1+a1*b3-a3*b1)*y+
((a2-b2)*c3+(b3-a3)*c2-a2*b3+a3*b2)*x+
((b1-a1)*c2+(a2-b2)*c1+a1*b2-a2*b1)*u3+
((a1-b1)*c3+(b3-a3)*c1-a1*b3+a3*b1)*u2+
((b2-a2)*c3+(a3-b3)*c2+a2*b3-a3*b2)*u1)/((b1-a1)*c2+(a2-b2)*c1+a1*
b2-a2*b1)
Give me a little while to verify this formula. Note this has a in it as well as u whereas the previous one did not. I may have posted the wrong formula originally.

[LATER]
That one looks right now. Not sure how the previous one got posted as it has no 'a' in it which means one whole point is missing or assumed to be a constant. Must have copied and pasted the wrong one. My Apologies.

If you want to verify here is the procedure i used.
First, using the three points A, B, C create two vectors that lie in the plane:
BA=B-A
CA=C-A

then calculate the cross product:
P=BA x CA

where P=[a,b,c].
This must be orthogonal to both vectors so it's a normal.

The equation of the plane is then:
a*i+b*j+c*k=0

Subst i=x-u1, j=y-u2, k=z-u3, then solve for z.
Note here a,b,c are the coefficients of i,j,k of P not the original a,b,c as points.

For reference, the cross product above is:
P=
((b1-a1)*(c2-a2)-(b2-a2)*(c1-a1))*k
+((b3-a3)*(c1-a1)-(b1-a1)*(c3-a3))*j
+((b2-a2)*(c3-a3)-(b3-a3)*(c2-a2))*i

and those are the points a,b,c.

Using the three points:
[[1,2,3],[3,2,1],[7,9,13]]

i get the same equation no matter which point i use for 'u':
z=(32*y-14*x-8)/14

I definitely copied the wrong formula from my notes before somehow, very sorry for that.
Substituting the points as before we get:
z=
(u1*((y2-y1)*z3+y1*z2+y3*(z1-z2)-y2*z1)
+x*((y1-y2)*z3+y3*(z2-z1)-y1*z2+y2*z1)
+y*((x2-x1)*z3+x1*z2+x3*(z1-z2)-x2*z1)
+u2*((x1-x2)*z3+x3*(z2-z1)-x1*z2+x2*z1)
+u3*((x2-x1)*y3+x1*y2+x3*(y1-y2)-x2*y1))
/((x2-x1)*y3+x1*y2+x3*(y1-y2)-x2*y1)

and that should work well too.

Last edited:
• AlbertHall

#### AlbertHall

Joined Jun 4, 2014
9,869
Hooray! Yes that one works (once I got it entered corectly!). It gives identical results to the other method and it doesn't matter which point I choose for u.
Now I just (!) have to convert them to program code. I think I will enter them both then I can check that their answers agree as a final check that I entered them correctly. Then I will choose one of them, probably on the basis of whichever is fastest.

#### MrAl

Joined Jun 17, 2014
7,591
Hooray! Yes that one works (once I got it entered corectly!). It gives identical results to the other method and it doesn't matter which point I choose for u.
Now I just (!) have to convert them to program code. I think I will enter them both then I can check that their answers agree as a final check that I entered them correctly. Then I will choose one of them, probably on the basis of whichever is fastest.

Hello again,

Oh that's great.
You know some of those formulas can be copy and pasted if they dont contain exponentiation because they are all in pure text. I dont know what language you are using though. Also, there may be a better factorization of those formulas as i did not look for every possibility.

I did one last check on the x1,y1,z1 form and was able to create the plane from three points as shown in this diagram. The plane came out as expected. #### AlbertHall

Joined Jun 4, 2014
9,869
The copy and paste is fine but all the values will be in arrays. The program will need to know where to find 'Z2', 'X1' etc.

#### MrAl

Joined Jun 17, 2014
7,591
The copy and paste is fine but all the values will be in arrays. The program will need to know where to find 'Z2', 'X1' etc.
Hi,

Oh ok. Maybe this will help then...

Here is a set of formulas in a form more suitable for programming. Note the partial symmetry.
A=(y2-y1)*z3+(y1-y3)*z2+(y3-y2)*z1
B=(x2-x1)*z3+(x1-x3)*z2+(x3-x2)*z1
C=(x2-x1)*y3+(x1-x3)*y2+(x3-x2)*y1
z=(u1-x)*A/C+(y-u2)*B/C+u3

and in terms of the differences this gets much more compact and easy to program:
A=(y21)*z3+(y13)*z2+(y32)*z1
B=(x21)*z3+(x13)*z2+(x32)*z1
C=(x21)*y3+(x13)*y2+(x32)*y1
z=(u1-x)*A/C+(y-u2)*B/C+u3

where we note the differences are sometimes repeated too.
The differences are:
y21=y2-y1
y13=y1-y3
x32=x3-x2
etc., so
xNM=xN-xM (and xM-xN would be WRONG)

That should make it a lot easier to program.
I checked the first set for equality with the previously posted long formula. The second set should be ok.

Note the somewhat non symmetry in u1-x and y-u2. I wanted to make them both the same but the signs worked out better in the constant B if i did it that way. Of course you could also switch to subtracting that term instead and 'fix' it back to u2-y instead.

The compact form for 'z' also shows that u3 appears as an offset which is interesting. For one thing if u3 was made the reference point as zero then it could be eliminated entirely, for what it is worth.

There's probably a matrix solution in there also but your program would have to have the ability to solve matrixes or at least handle them in some way.

Last edited: