fast geometric proximity predicate

Question

I have 3 points (A, B and X) and a distance (d). I need to make a function that tests if point X is closer than distance d to any point on the line segment AB.

The question is firstly, is my solution correct and then to come up with a better (faster) solution.

My first pass is as follows

AX = X-A
BX = X-B
AB = A-B

    // closer than d to A (done squared to avoid needing to compute the sqrt in mag)
If d^2 > AX.mag^2  return true

    // closer than d to B
If d^2 > BX.mag^2 return true

    // "beyond"  B
If (dot(BX,AB) < 0) return false

    // "beyond"  A
If (dot(AX,AB) > 0) return false

    // find component of BX perpendicular to AB
Return (BX.mag)^2 - (dot(AB,BX)/AB.mag)^2 < d^2

This code will end up being run for a large set of P's and a large set of A/B/d triplets with the intent of finding all P's that pass for at least one A/B/d so I suspect that there is a way to reduce overall the cost based on that but I haven't looked into that yet.

(BTW: I am aware that some reordering, some temporary values and some algebraic identities could decrease the cost of the above. I just omitted them for clarity.)

---

EDIT: this is a 2D problem (but solution that generalizes to 3D would be cool

Edit: On further reflection, I expect the hit rate to be around 50%. The X point can be formed in a nested hierarchy so I expect to be able to prune large subtrees as all-pass and all-fail. The A/B/d limiting the triplets will be more of a trick.

Edit: d is in the same order of magnitude as AB.

---

edit: Artelius posted a nice solution. I'm not sure I understand exactly what he's getting at as I got off on a tangent before I fully understood it. Anyway another thought came to mind as a result:

• First Artelius' bit, pre-cacluate a matrix that will place AB centered ate the origin and aligned with the X-axis. (2 adds, 4 muls and 2 adds)
• fold it all into the 1st quadrant (2 abs)
• scale in X&Y to make the central portion of the zone into a unit square (2 mul)
• test if the point is in that square (2 test) is so quit
• test the end cap (go back to the unscaled values
  • translate in x to place the end at the origin (1 add)
  • square and add (2 mul, 1 add)
  • compare to d^2 (1 cmp)

I'm fairly sure this beats my solution.

(if nothing better comes along sone Artelius gets the "prize" :)

Answer 1

Hmmmmmmm.... What's the hit-rate? How often does point "X" meet the proximity requirements?

I think your existing algorithm is good, and any additional optimizations will come from tuning to the real data. For example, if the "X" point meets the proximity test 99% of the time, then I think your optimization strategy should be considerably different than if it only passes the test 1% of the time.

---

Incidentally, when you get to the point where you're running this algorithm with thousands of points, you should organize all the points into a K-Dimensional Tree (or KDTree). It makes the calculation of "nearest-neighbor" much simpler.

Your problem is a little more complex than a basic nearest-neighbor (because you're checking proximity-to-a-line-segment rather than just proximity-to-a-point) but I still think the KDTree will be handy.

Answer 2

This code will end up being run for a large set of P's and a large set of A/B/d triplets with the intent of finding all P's that pass for at least one A/B/d so I suspect that there is a way to reduce overall the cost based on that but I haven't looked into that yet.

In the case when d ~ AB, for a given point X, you can first test whether X belongs to one of the many spheres of radius d and center Ai or Bi. Look at the picture:

     ......        .....
  ...........   ...........
 ...........................
.......A-------------B.......
 ...........................
  ...........   ...........
     .....         .....

The first two tests

If d^2 > AX.mag^2 return true
If d^2 > BX.mag^2 return true

are the fastest ones, and if d ~ AB they are also the ones with highest probability to succeed (given that the test succeeds at all). Given X, you can do all the "sphere tests" first, and then all the final ones:

Return (BX.mag)^2 - (dot(AB,BX)/AB.mag)^2 

Answer 3

If I read this correctly, then this is almost the same as the classic ray/sphere intersection test as used in 3D ray-tracing.

In this case you have a sphere at location X of radius d, and you're trying to find out whether the line AB intersects the sphere. The one difference with ray-tracing is that in this case you've got a specific line AB, whereas in ray-tracing the ray is normally generalised as as origin + distance * direction, and you don't care how far along the infinite line AB+ it is.

In pseudo-code from my own ray-tracer (based on the algorithm given in "An Introduction to Ray-tracing" (ed. Glassner):

Vector v = X - A
Vector d = normalise(B - A)  // unit direction vector of AB
double b = dot(v, B - A)
double discrim = b^2 - dot(v, v) + d^2
if (discrim < 0)
     return false            // definitely no intersection

If you've got this far, then there's some chance that your condition is met. You just have to figure out whether the intersection(s) is on the line AB:

discrim = sqrt(discrim)
double t2 = b + discrim
if (t2 <= 0)
    return false             // intersection is before A

double t1 = b  - discrim

result = (t1 < length(AB) || (t2 < length(AB))

Answer 4

If your set of (A,B,d) in fixed, you can calculate a pair of matrices for each to translate the co-ordinate system, so that the line AB becomes the X axis, and the midpoint of AB is the origin.

I think this is a simple way to construct the matrices:

trans = - ((A + B) / 2)        // translate midpoint of AB to origin

rot.col1 = AB / AB.mag         // unit vector in AB direction

                        0 -1    
rot.col2 = rot.col1 * (      ) // unit vector perp to AB
                        1  0 

rot = rot.inverse()            // but it needs to be done in reverse

Then you just take each X and do rot * (X + trans).

The region in question is actually symmetric in both the x and y axes now, so you can take the absolute value of the x co-ordinate, and of the y co-ordinate.

Then you just need to check:

y < d && x < AB.mag/2            //"along" the line segment
|| (x - AB.mag/2)^2 + y^2 < d^2  // the "end cap".

You can do another trick; the matrix can scale down by a factor of AB.mag/2 (remember the matrices are only calculated once per (A,B) meaning that it's better that finding them is slower, than the actual check itself). This means your check becomes:

y < 2*d/AB.mag && x < 1
|| (x - 1)^2 + y^2 < (2*d/AB.mag)^2

Having replaced two instances of AB.mag/2 with the constant 1, it might be a touch faster. And of course you can precalculate 2*d/AB.mag and (2*d/AB.mag)^2 as well.

Whether this will work out faster than other methods depends on the inputs you give it. But if the length of AB is long compared to d, I think it will turn out considerably faster than the method you posted.

Answer 5

If the # of sets of A/B/d are large, and you're definitely in 2D, consider using R-trees (or their octagonal equivalents) where each entry in the R-tree is the minimum bounding box of the A/B/d triple. This would let you eliminate having to test a lot of the A/B/d triples & focus your CPU power only on the few ones where each point X is within the bounding boxes of the A/B/d triple. Then do a more detailed test like the one you mention.