6
$\begingroup$

I need to compute the distance between two line segments in a project. After googling, I found this algorithm and used it to implement a Mathematica version:

distSegToSeg[{p1_, p2_}, {q1_, q2_}] := Module[{small = 10^(-6), u, v, w, a, b, c, d, e, D, sc, sN, sD, tc, tN, tD}, u = p2 - p1; v = q2 - q1; w = p1 - q1; a = u.u; b = u.v; c = v.v; d = u.w; e = v.w; D = a*c - b*b; sD = D; tD = D; If[D < small, sN = 0; sD = 1; tN = e; tD = c, sN = b*e - c*d; tN = (a*e - b*d); If[sN < 0, sN = 0.0; tN = e; tD = c, If[sN > sD, sN = sD; tN = e + b; tD = c;]]]; If[tN < 0, tN = 0; If[-d < 0, sN = 0, If[-d > a, sN = sD, sN = -d; sD = a]], If[tN > tD, tN = tD; If[-d + b < 0, sN = 0, If[-d + b > a, sN = sD, sN = -d + b; sD = a]]]]; sc = If[Norm[sN] < small, 0, sN/sD]; tc = If[Norm[tN] < small, 0, tN/tD]; N[Norm[w + sc*u - tc*v]] ]

The code is not in good Mathematica style, and it's very slow. Do you have any idea to make this faster? Do you have a better way to compute the distance between two line segments in 3D?

Improve this question
$\endgroup$
0

5 Answers 5

Reset to default
5
$\begingroup$

You can easily see that if the nearest approach of the two infinite lines does not fall within both segments then you must calculate the nearest distance of all four end points to the other segment.

 pointsegdis[{seg_, pointlist_}] := Module[{u = Subtract @@ seg, mean = Plus @@ seg/2}, {#, (mean - u Sign@# Min[1, Abs@#]/2) &@(-(2 ( # - mean).u )/u.u) } & /@ pointlist ] segsegdis[s1_, s2_] := Module[{u, v, w, d, t1, t2}, {u, v, w} = Subtract @@ # & /@ {Sequence @@ #, First /@ #} &@{s1, s2}; If[(d = u.u v.v - (v.u)^2) != 0 && Abs[t2 = 2 (u.u v.w - u.w v.u)/d + 1] <= 1 && Abs[t1 = 2 (v.u v.w - v.v u.w)/d + 1] <= 1, Plus @@ (#[[1]] (1 + {1, -1} #[[2]] )) /2 & /@ {{s1, t1}, {s2, t2}} , First@SortBy[Flatten[pointsegdis /@ {{s1, s2}, {s2, s1}}, 1], Total@((Subtract @@ #)^2) &] ] ] ; {s1, s2} = RandomReal[{-1, 1}, {2, 2, 3}]; connect = segsegdis[s1, s2]; Graphics3D[{Line /@ {s1, s2}, {Red, Line[connect]}}]

enter image description here

This is several orders of magnitude faster than Minimise or FindMinimum

 testcases = RandomReal[{-1, 1}, {100, 2, 2, 3}]; (aa = Norm[ Subtract @@ segsegdis[#[[1]], #[[2]]] ] & /@ testcases) //Timing // First (bb = First@ Minimize[{ Norm[(#[[1, 1]] (1 - t1)/2 + #[[1, 2]] (1 + t1)/2) - (#[[2, 1]] (1 - t2)/2 + #[[2, 2]] (1 + t2)/2) ] , {-1<=t1 <= 1 , -1 <= t2 <= 1}}, {t1, t2}] & /@ testcases) // Timing //First (cc = First@ FindMinimum[{ Norm[(#[[1, 1]] (1 - t1)/2 + #[[1, 2]] (1 + t1)/2) - (#[[2, 1]] (1 - t2)/2 + #[[2, 2]] (1 + t2)/2) ] , {-1 <= t1 <= 1 , -1 <= t2 <= 1}}, {t1, t2}] & /@ testcases) // Timing // First

{0.015600, 12.620481, 2.714417}

verify all three give the same result to reasonable numerical precision

 Max[Abs[aa - bb]] -> 5.82532*10^-10 Max[Abs[aa - cc]] -> 2.25578*10^-6
Improve this answer
$\endgroup$
2
$\begingroup$

The link you gave shows how to find the distance between any two points on the lines. You can then use Mathematica's Minimize function to find the shortest distance. Like this.

minDist[{p1_, p2_}, {q1_, q2_}] := Module[{P, Q, u, v, w}, u = Normalize[p2 - p1]; v = Normalize[q2 - q1]; P[s_] := p1 + s u; Q[t_] := q1 + t v; w[s_, t_] := P[s] - Q[t]; First[Minimize[Norm[w[s, t]], {s, t}]] ]

Here, u and v are unit vectors pointing from p1 to p2 and q1 to q2, respectively. The functions P and Q represent the two lines. The function w is a vector going between P[s] and Q[t]. Minimize the length of that vector.

Improve this answer
$\endgroup$
1
  • $\begingroup$ It's much clear, but I test it , it's much slower than my previous code $\endgroup$ Commented Apr 4, 2014 at 0:32
2
$\begingroup$

If you want a method that uses built-in functions, this was suggested to me by a colleague 

distLineToLine[line1_, line2_] := Module[{distFunc}, distFunc = RegionDistance[Line@line1]; NMinValue[distFunc[{x, y, z}], {x, y, z} ∈ Line@line2] ]

But it isn't so fast,

lines = {{{24.134, -147.93, 193.37}, {-9.1854, -85.421, 96.284}}, {{-135.82, -104.67, 38.896}, {-177.11, -32.589, -74.348}}}; distLineToLine @@ lines // RepeatedTiming (* {0.17, 140.358} *)

The best answer here, in my opinion, is the code from the OP, which is about an order of magnitude faster than the accepted answer, although the code doesn't look very pretty. You can speed it up even more by compiling it (code pasted to gist to save space here)

distSegToSegComp = << "https://git.io/vDJMK"

Now compare the timing,

distSegToSeg @@ lines // RepeatedTiming segsegdis @@ lines // Apply[Subtract] // Norm // RepeatedTiming distSegToSegComp @@ lines // RepeatedTiming (* {0.000054, 140.358} *) (* {0.00012, 140.358} *) (* {3.7*10^-6, 140.358} *)

So the compiled version wins hands down.

Improve this answer
$\endgroup$
1
$\begingroup$

Another way of doing this (Link) is to find the mutual perpendicular between the two lines using the cross product, converting this to a unit vector, and then using the dot product between that cross product, and any vector going between the two lines. Like this:

newMinDist[{p1_, p2_}, {q1_, q2_}] := Module[ {u, v, n, w}, u = p2 - p1; v = q2 - q1; n = Normalize[Cross[u,v]]; w = q1 - p1; Dot[w,n] ]

Again, u and v are vectors headed along the lines. The vector n is a unit vector in the direction of the cross product of u and v (the unit vector normal to both u and v). The w is any vector between the ps and the qs. You could swap any value of 1 or 2 here (for instance q2-p1). Then w.n gives the shortest distance between the lines.

Gives the same result, and Timing shows this method is better than an order of magnitude faster than the previous one.

Improve this answer
$\endgroup$
0
$\begingroup$

Suppose the equation of two lines are

(x-x1)/A1=(y-y1)/B1=(z-z1)/C1 (x-x2)/A2=(y-y2)/B2=(z-z2)/C2

Then the oriented distance between them is

Det[{{x2-x1,y2-y1,z2-z1},{A1,B1,C1},{A2,B2,C2}}]/Norm[Cross[{A1,B1,C1},{A2,B2,C2}]]

You may need Abs to get the distance

Improve this answer
$\endgroup$
2
  • $\begingroup$ what I want is the distance of two line segment, not two infinite line $\endgroup$ Commented Apr 4, 2014 at 4:00
  • $\begingroup$ @tintin its the same calculation now just check that your on the segment or the endpoint is closest. $\endgroup$ Commented Apr 5, 2014 at 6:54

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.