5

While looking at various methods for point-in-triangle testing (2D case), I found that the method which uses barycentric coordinates is the most used one. Here is a StackOverflow answer which explains it.

Why is this method the most preferred one? It probably has to do with doing less calculations, but what about numerical stability? Is this algorithm better suited than say, the "same side" technique, for cases in which the point is particularly near the border?

Community
  • 1
  • 1
rubik
  • 8,814
  • 9
  • 58
  • 88
  • I'm not sure about robustness. But in many cases you need the barycentric coordinates for the next step, anyway. E.g. in Raytracing, you check if a ray hits the triangle. And if it hits, interpolate colors from the triangle's vertices using the barycentric coordinates. – Nico Schertler May 25 '16 at 06:10

1 Answers1

3

If you solve it:

p = p0 + (p1 - p0) * s + (p2 - p0) * t
s = <0.0,1.0>
t = <0.0,1.0>
s+t<=1.0

While solwing this system:

p.a = p0.a + (p1.a - p0.a) * s + (p2.a - p0.a) * t
p.b = p0.b + (p1.b - p0.b) * s + (p2.b - p0.b) * t
----------------------------------------------------

You got two algebraic options:

I.  t = (p.a - p0.a - (p1.a - p0.a) * s) / (p2.a - p0.a)
II. p.b = p0.b + (p1.b - p0.b) * s + (p2.b - p0.b) * t
----------------------------------------------------
II. p.b = p0.b + (p1.b - p0.b) * s + (p2.b - p0.b) * (p.a - p0.a - (p1.a - p0.a) * s) / (p2.a - p0.a)
II. s = (p.b-p0.b) / ( (p1.b-p0.b) + ( (p2.b-p0.b)*(p.a-p0.a-(p1.a-p0.a)/(p2.a-p0.a) ) )
...

And:

I.  s = (p.a - p0.a - (p2.a - p0.a) * t) / (p1.a - p0.a)
II. p.b = p0.b + (p1.b - p0.b) * s + (p2.b - p0.b) * t
----------------------------------------------------
...

Which gives you 2 algebraic solution options. To ensure stability you should divide with big enough magnitudes. So you should choose the axises (a,b -> x,y) , and point order so you are not dividing by zero or small magnitude numbers.

To avoid this you can use matrix approach

p.a = p0.a + (p1.a - p0.a) * s + (p2.a - p0.a) * t
p.b = p0.b + (p1.b - p0.b) * s + (p2.b - p0.b) * t
--------------------------------------------------
|p.a|   | (p1.a - p0.a) , (p2.a - p0.a) , p0.a |   | s |
|p.b| = | (p1.b - p0.b) , (p2.b - p0.b) , p0.b | * | t |
| 1 |   |       0       ,       0     ,     1  |   | 1 |
--------------------------------------------------------
| s |           | (p1.a - p0.a) , (p2.a - p0.a) , p0.a |   | p.a |
| t | = inverse | (p1.b - p0.b) , (p2.b - p0.b) , p0.b | * | p.b |
| 1 |           |       0       ,       0       ,   1  |   |  1  |
------------------------------------------------------------------

Also here you got more options for axises order, point order so that the inverse matrix is computable. If you use sub-determinant approach for inverse matrix solution then the only thing that should matter is the final division step. So you can choose the orders until you have nonzero determinant.

Spektre
  • 49,595
  • 11
  • 110
  • 380
  • I'm diving in the numerical stability topic, I just have one more question. You say that to preserve stability we should divide by big enough magnitudes. Shouldn't we worry only about close-to-zero denominators? For example, if the denominators are 4, 2, 5, there would be no difference between them as they're big enough. Is that correct? – rubik May 28 '16 at 06:58
  • @rubik I am no expert in the matter but from my experience it is always better to use as big divider as you can Otherwise you can get oscillations around bounds (in critical conditions) making pixelated border artifacts. But that is only in rare circumstances. If you are on floating point variables then you should be fine. – Spektre May 28 '16 at 08:05