Playing with algorithms.

[MHD]

I like thinking of algorithms and how to use the fastest atomic operations and the fewest and so on...
The other day it occurred to me an intuitive way to test if a point is inside a triangle.
(Assume all numbers floating points of some kind)

The euclidean distance is relatively expensive, because it uses square root. We can alleviate that by squaring all our distances instead.

Code: Select all

let dist(P1, P2) = x * x + y * y
    where x = P1.x - P2.x
          y = P1.y - P2.y

Three additions, two multiplications.

When we look at a triangle, all the points inside the triangle have distances to the corners shorter than or equal to the sides.

Code: Select all

let inside(P, A, B, C) = pa < ab && pa < ca && pb < ab && pb < bc && pc < bc && pc < ca
    where pa = dist(P, A)
          pb = dist(P, B)
          pc = dist(P, C)
          ab = dist(A, B)
          bc = dist(B, C)
          ca = dist(C, A)

Six floating comparisons, five boolean conjunctions, six uses of dist (18 additions, 12 multiplications).

This algorithm just occurred to me. I know this is probably reinventing the wheel, but hey, it's my hobby.
Most of the approaches I see are to use vector operations which I can see makes sense in 3d and when you have a GPU handy. What if you don't?

Do you ever discover or "play" with algorithms?
What novel approaches have you encountered?
Got any gems to share?

[phlip]

I'm pretty sure that doesn't work... consider the triangle (0,0), (1,0), (0,1) and the point (0.6, 0.6). That satisfies your criterion, but is outside the triangle. Meanwhile, consider the triangle (0,0), (1,0), (0,2) and the point (0.02, 1.9). That point is outside the triangle, but fails your check.

[MHD]

I'm pretty sure that doesn't work... consider the triangle (0,0), (1,0), (0,1) and the point (0.6, 0.6). That satisfies your criterion, but is outside the triangle. Meanwhile, consider the triangle (0,0), (1,0), (0,2) and the point (0.02, 1.9). That point is outside the triangle, but fails your check.

This should fix it...

Code: Select all

let inside(P, A, B, C) = pab > ab && pab < (bc + ca) && pac > ac && pac < (ab + bc) && pbc > bc && pbc < (ab + ca)
    where pa = dist(P, A)
          pb = dist(P, B)
          pc = dist(P, C)
          ab = dist(A, B)
          bc = dist(B, C)
          ca = dist(C, A)
          pab = pa + pb
          pac = pa + pc
          pbc = pb + pc

[mfb]

No, triangles do not work like that.

pab>ab is true iif P is not on the edge AB.
pab < bc+ca is true within an ellipse with A and B as focal points and C on its line.

So phlips first triangle and point still get the wrong result with your method.

A strong hint: Instead of checking distances, it is better to check on which side of a line your point P is.
In the worst case, you need 15 subtractions, 6 multiplications and 3 conditional jumps .

[Yakk]

Naive:
Norm (B-A) . (p-A) > 0
Norm (C-B) . (p-B) > 0
Norm (A-C) . (p-C) > 0
where Norm returns a vector rotated 90 degrees in the appropriate direction, and ABC is oriented in the appropriate direction (if the orientations aren't correct, you end up with an anti-triangle, heh), and . is the dot product operator.

That is 15 additions, 6 multiplications, and 3 if clauses, so is probably mfb's solution.

Norm is surprisingly efficient, requiring a sign shift and a swap.

I can strip 2 subtractions out of there by using the same base point for two of the sides.

As another advantage, this can be generalized to a convex n-gon (or a simplex, which is always convex)

[mr-mitch]

In a 2D problem, if you use the cross product instead, you don't need to take the normal vector. There is also an equal number of operations as in the dot product.

(a,b) x (c,d) = ad - bc
(a,b) . (c,d) = ac + bd
Although, of course, taking the norm and then taking the dot product gives you the same thing.
Basically you can avoid the shift sign and swap.

Instead the operations you listed all need to be the same sign (pointing in the same direction). If you store your simple polygon in clockwise order, then each product will need to be positive (or non-negative if on the line counts).

(A_i+1 - A_i) x (p - A_i) > 0

[korona]

In a 2D problem, if you use the cross product instead, you don't need to take the normal vector. There is also an equal number of operations as in the dot product.

(a,b) x (c,d) = ad - bc
(a,b) . (c,d) = ac + bd
Although, of course, taking the norm and then taking the dot product gives you the same thing.
Basically you can avoid the shift sign and swap.

What?

||(1, 2)|| * ||(3, 4)|| = sqrt(5) + 5 (taking the norm and then the dot product? (the "dot product" on R is simple multiplication?))
(1, 2) x (3, 4) = 4 - 6 = -2 (i've never seen this operation called cross product)
(1, 2) . (3, 4) = 3 + 6 = 9

These operations are not related in any way (??). Computing normals in 2D is easy as (-b, a) is a normal of (a, b).
Are you sure you meant what you wrote?

[Yakk]

Let o() be the (well, a) orthogonal-generating function (ie, o(x,y) -> (y, -x)).

Examine
a x b
and
a x o(b)

In particular:
(1, 2) . o(3, 4) = (1,2) . (4, -3) = 4-6 = -2

The other one (in R^2) simply does a sign-reversal, and corresponds to the other 90 degree rotation (and the other way of embedding the cross-product result back into the original space in R^2, but that is getting rather obscure).

[mr-mitch]

Yes, it's not really the cross product. In fact, in a 2D problem, the operation is in fact ((x,y,0) x (a,b,0)) . k. By the norm I was referring to Yakk's post, finding the orthogonal vector (normal to).

You can generalise this to a planar surface by making k your normal vector to the plane.
In 3D problems, you would extend it to making sure it's on the same side of each face/plane.

Basically what you're doing is you're making sure each vector from the cross product points in the same direction. In a 2D problem it's better than using the dot product because you don't have to find the normal of the vector, not that it is hard to do anyway. It gives you the same result.

Norm(a,b) . (c,d) = (-b, a) . (c,d) = ad - bc
(a,b,0) x (c,d,0) . k = (0,0, ad - bc) . k = ad - bc.

It depends on how define Norm. If they all have the same sign then it's inside (all above or below the surface). If you allow for 'on the line', then zeroes are wild cards.