# Fast comparison of line segments lengths

I have two line segments given by their endpoints $(a_1,a_2)$, $(b_1,b_2)$ in $R^3$ and want to know if they have the same length (up to some error), so that the naive test looks like

$$|\, \Vert a_1-a_2\Vert- \Vert b_1-b_2\Vert\, | < d_\text{threshold}$$

However, I need to test a large set of segments so I would like to get rid of the sqrt within the norm computation. Each segment only shows up once in the full computation so that I can precompute and store lenghts.

If $a$ is the squared norm of $a_1-a_2$, I can write the equation as

$$|\sqrt{a} - \sqrt{b}| < d_\text{threshold}$$

If I square both sided, I only get a

$$a^2 - 2 \sqrt{a} \sqrt{b} + b^2 < d_\text{threshold}^2$$

which again forces me to compute the roots.

Am I missing something obvious?

You could also rewrite the expression to avoid square roots entirely. Reorder $a,b$ so that $a>b>0$, and let the inequality be $$\sqrt{a}-\sqrt{b}<\tau.$$ This is equivalent to $$a < \tau^2 + 2\tau\sqrt{b} + b \qquad\Leftrightarrow\qquad 2\tau\sqrt{b} > a-\tau^2-b,$$ which is in turn equivalent to $$\Big( a < \tau^2+b \quad\text{or}\quad 4\tau^2b>(a-\tau^2-b)^2 \Big)$$ which does not involve square roots.

• That looks really good. I just don't understand how you got the (a<t**2+b). Could you clarify it a bit before I accept the answer? – FooBar Feb 15 '16 at 19:17
• @FooBar If $x=2\tau\sqrt{b}$ is known to be positive then $x>y$ holds if and only if $y$ is negative or $y$ is positive and $x^2>y^2$. Otherwise the two lines wouldn't be equivalent. – Kirill Feb 15 '16 at 19:28

Generally speaking, the $\mathtt{sqrt}$ function is going to be the slowest part of that. Fortunately, if their lengths are the same, then the squares of their lengths are also the same. Therefore, you can just omit that part and compare the squares of the norm.

I don't know of any way to make that go too much faster using cross-platform code.

However, if you really want to get it running quickly, you can always take a look at what processor you're running on and find out if it supports vector instructions (ie, SSE or AVX) and use those intrinsics. https://software.intel.com/sites/landingpage/IntrinsicsGuide/

By doing this, you can take advantage of instruction-level parallelism in the computation and achieve huge speedups, albeit at the cost of your time. It's possible that your compiler is already doing this for you, however, especially if you compiler with (for gcc) -march=native. You should check the assembly that comes out of the compiler before spending time vectorizing the code yourself.

• The squares of the lengths are similar, but how does the test for the squared lengths look like? If I compare the absolute difference of the squared lengths, then longer segments are less probable to be accepted. – FooBar Feb 15 '16 at 17:27
• If you have an acceptance condition that is related to the absolute differences of the lengths, then you may not be able to avoid the cost of computing these lengths. I think you'll have to determine if the difference of squares is sufficient for your application. – Tyler Olsen Feb 15 '16 at 18:29

I want to submit a separate answer to this question, because there is a way to automatically simplify such formulas to remove square roots. I will use sage and QEPCAD (both free). The tutorial is here.

The inequality $$|\sqrt{a}-\sqrt{b}|<t$$ is equivalent to the statement that there exist $u, v$ such that $$\exists u,v: \qquad (u-v)^2 < t^2, \quad u^2 = a, \quad\text{and}\quad v^2 = b.$$ These are quantified polynomial in-/equalities and can be simplified with QEPCAD using the following simple session in sage:

var('a,b,u,v,t')
qepcad(qf.E([u,v], qf.and_((u-v)^2 < t^2, u^2 == a, v^2 == b)))


The output gives the equivalent conditions for the original inequality to hold:

a >= 0 /\ b >= 0 /\ [ t^4 - 2 b t^2 - 2 a t^2 + b^2 - 2 a b + a^2 < 0 \/ t^2 - b - a > 0 ]


In other words, it is equivalent to $$(t^4 - 2 b t^2 - 2 a t^2 + b^2 - 2 a b + a^2 < 0) \lor (t^2 - b - a > 0).$$

Instead of putting a threshold on the differences, put a threshold on the ratio:

$$\frac{ \lVert \mathbf{b}_2-\mathbf{b}_1 \rVert}{ \lVert \mathbf{a}_2-\mathbf{a}_1 \rVert}< \tau$$ You could then reduce this to: $$r =\frac{ \lVert \mathbf{b}_2-\mathbf{b}_1 \rVert^2}{ \lVert \mathbf{a}_2-\mathbf{a}_1 \rVert^2}< \tau^2$$

and get rid of the square root. If you are disturbed by the fact that $r$ depends on the order selected, than you could consider a short if statement.

• I also thought about this way around the problem, but this is not a solution to my problem. My error does not depend on the distance, but is absolute. – FooBar Feb 15 '16 at 17:24
• You could always scale your $\tau$ with the absolute distance you have, to convert it into this form. – Tolga Birdal Feb 16 '16 at 9:18

If $a$ is the squared norm of $a_1-a_2$, then the equation can written as:

$$|\sqrt{a} - \sqrt{b}| < d_\text{threshold}$$

But squaring both sides results in:

$$a - 2 \sqrt{a} \sqrt{b} + b < d_\text{threshold}^2$$ or $$- 2 \sqrt{a} \sqrt{b} < d_\text{threshold}^2 - (a + b)$$

so $$(a + b - d_\text{threshold}^2)< 2\sqrt{a} \sqrt{b}$$

Let $c = a + b - d_\text{threshold}^2$

Thus $$c < 2\sqrt{a} \sqrt{b}$$

If c >= 0 then the inequality becomes $$c^2 < 4ab$$ Otherwise use $$c^2 > 4ab$$

The inequality flips due to effectively multiplying by the negative value $c$. From a programming standpoint, I like the fact that very few operations would be required to determine inclusion.

The original answer neglected the fact that $d^2$ could be larger than $a + b$.

• But with this approach, I think you would lose the metric intuition. Because the threshold is not on the difference of the square roots and not actually the distances themselves. Am I right? – Tolga Birdal Feb 16 '16 at 9:19
• It really depends on what the goals are. Sometimes intuition must give way to computational speed. It helps if you are able to convince yourself that the equation/inequality is grounded in reality. – Robert Benson Feb 16 '16 at 13:48
• This seems to incorrectly reject $(a,b,d) = (\frac1{16},\frac1{32},\frac12)$. – Kirill Feb 23 '16 at 0:09
• @Kirill I see the problem. It comes from the fifth line. The LHS is always < 0, but the RHS may or may not be depending on the particular inputs. I think I see a workaround that would just need to do one more comparison to tell how the inequality should go. – Robert Benson Feb 23 '16 at 20:04
• @RobertBenson Yes, compare with: scicomp.stackexchange.com/a/23179/713 – Kirill Feb 23 '16 at 20:17