## Playing with algorithms.

A place to discuss the science of computers and programs, from algorithms to computability.

Formal proofs preferred.

Moderators: phlip, Moderators General, Prelates

### Playing with algorithms.

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`

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?
EvanED wrote:be aware that when most people say "regular expression" they really mean "something that is almost, but not quite, entirely unlike a regular expression"

MHD

Posts: 631
Joined: Fri Mar 20, 2009 8:21 pm UTC
Location: Denmark

### Re: Playing with algorithms.

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.
While no one overhear you quickly tell me not cow cow.

phlip
Restorer of Worlds

Posts: 6964
Joined: Sat Sep 23, 2006 3:56 am UTC
Location: Australia

### Re: Playing with algorithms.

phlip wrote: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`
EvanED wrote:be aware that when most people say "regular expression" they really mean "something that is almost, but not quite, entirely unlike a regular expression"

MHD

Posts: 631
Joined: Fri Mar 20, 2009 8:21 pm UTC
Location: Denmark

### Re: Playing with algorithms.

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 .
mfb

Posts: 824
Joined: Thu Jan 08, 2009 7:48 pm UTC

### Re: Playing with algorithms.

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)
One of the painful things about our time is that those who feel certainty are stupid, and those with any imagination and understanding are filled with doubt and indecision - BR

Last edited by JHVH on Fri Oct 23, 4004 BCE 6:17 pm, edited 6 times in total.

Yakk
Poster with most posts but no title.

Posts: 10212
Joined: Sat Jan 27, 2007 7:27 pm UTC
Location: E pur si muove

### Re: Playing with algorithms.

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).

(Ai+1 - Ai) x (p - Ai) > 0
mr-mitch

Posts: 473
Joined: Sun Jul 05, 2009 6:56 pm UTC

### Re: Playing with algorithms.

mr-mitch wrote: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?
korona

Posts: 241
Joined: Sun Jul 04, 2010 8:40 pm UTC

### Re: Playing with algorithms.

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).
One of the painful things about our time is that those who feel certainty are stupid, and those with any imagination and understanding are filled with doubt and indecision - BR

Last edited by JHVH on Fri Oct 23, 4004 BCE 6:17 pm, edited 6 times in total.

Yakk
Poster with most posts but no title.

Posts: 10212
Joined: Sat Jan 27, 2007 7:27 pm UTC
Location: E pur si muove

### Re: Playing with algorithms.

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.
mr-mitch

Posts: 473
Joined: Sun Jul 05, 2009 6:56 pm UTC