Re: Distance point <=> straight line in space

Daniel Kraft <d@xxxxxxxx> writes:

in a program, I've got a straight line in space defined by two points
(p1 and p2) on it, and a set (~ 10000) points q_i in space. Then I
need to calculate for each of those points the distance to the
straight. At the moment I'm doing it like this:

vector g = normalized (p2 - p1) // Find normalized direction
for i = 1 to n
distance_i = norm (cross_product (p2 - q_i, g))

So this costs (I think):

* +- sqrt
3 [from p2 - q_i]
6 5 [cross product]
3 2 1 [norm]

or 9 multiplies, 10 add/subtracts and 1 square root.

I believe this formula is fairly straight-forward (a single
cross-product and a norm); unfortunatelly, my program spents a large
amount of time (~ 95% of total runtime) in this routine and the
current duration is unacceptable.

Is there a better approach I could take to do this calculation that
requires less operations?

The uglier formula[1] involves pre-computing

V = (p2 - p1)/norm_squared(p2 - p1);

where norm_squared (x, y, z) is x^2 + y^2 + z^2. Then, for each
point, the distance is:

sqrt(norm_squared(Q) - (Q.V)^2) where Q = p1 - q_i;

* +- sqrt
3 [from p1 - q_i]
3 2 [norm_squared(Q)]
3 2 [dot product]
1 1 [square and subtract]
1 [sqrt]

which I make to be 7 multiplies, 8 add/subtracts and a square root.
Of course, I may have messed up the algebra.

... The output of gcc -O3 seems already
reasonably optimized, i.e., the cross-product and norm calculation is
fully inlined. I don't think hand-optimized inline-assembly could
help here much, right?

If there is any reason to transform your co-ordinate system, you might
want to shift and rotate the world so that p1 and p2 are on the z-axis
(say). Then all you distances are 2D from the origin. This feels
more expensive unless there is some other payoff to changing the
co-ordinate system.

[1] See, for example:
http://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html

--
Ben.
.