Stephen Royle

An occasional series in esoteric programming issues.
As part of a larger analysis project I needed to implement a short program to determine the closest distance of two line segments in 3D space. This will be used to sort out which segments to compare&#8230; like I say, part of a bigger project. The best method to do this is to find the closest distance one segment is to the other when the other one is represented as an infinite line. You can then check if that point is beyond the segment if it is you use the limits of the segment to calculate the distance. There&#8217;s a discussion on stackoverflow <a href="http://stackoverflow.com/questions/627563/calculating-the-shortest-distance-between-two-lines-line-segments-in-3d">here</a>. The solutions point to <a href="http://geomalgorithms.com/a07-_distance.html">one in C++</a> and one in <a href="http://stackoverflow.com/a/702174">MATLAB</a>. The C++ version is easiest to port to Igor due to the similarity of languages, but the explanation of the MATLAB code was more approachable. So I ported that to Igor to figure out how it works.
The MATLAB version is:
<pre class="brush: matlabkey; title: ; notranslate" title="">
&gt;&gt; P = &#x5B;-0.43256 -1.6656 0.12533]
P =
 -0.4326 -1.6656 0.1253
&gt;&gt; Q = &#x5B;0.28768 -1.1465 1.1909]
Q =
 0.2877 -1.1465 1.1909
&gt;&gt; R = &#x5B;1.1892 -0.037633 0.32729]
R =
 1.1892 -0.0376 0.3273
&gt;&gt; S = &#x5B;0.17464 -0.18671 0.72579]
S =
 0.1746 -0.1867 0.7258
&gt;&gt; N = null(P-Q)
N =
 -0.3743 -0.7683
 0.9078 -0.1893
 -0.1893 0.6115
&gt;&gt; r = (R-P)*N
r =
 0.8327 -1.4306
&gt;&gt; s = (S-P)*N
s =
 1.0016 -0.3792
&gt;&gt; n = (s - r)*&#x5B;0 -1;1 0];
&gt;&gt; n = n/norm(n);
&gt;&gt; n
n =
 0.9873 -0.1587
&gt;&gt; d = dot(n,r)
d =
 1.0491
&gt;&gt; d = dot(n,s)
d =
 1.0491
&gt;&gt; v = dot(s-r,d*n-r)/dot(s-r,s-r)
v =
 1.2024
&gt;&gt; u = (Q-P)'\((S - (S*N)*N') - P)'
u =
 0.9590
&gt;&gt; P + u*(Q-P)
ans =
 0.2582 -1.1678 1.1473
&gt;&gt; norm(P + u*(Q-P) - S)
ans =
 1.0710
</pre>
and in IgorPro:
<pre class="brush: plain; title: ; notranslate" title="">
Function MakeVectors()
	Make/O/D/N=(1,3) matP={{-0.43256},{-1.6656},{0.12533}}
	Make/O/D/N=(1,3) matQ={{0.28768},{-1.1465},{1.1909}}
	Make/O/D/N=(1,3) matR={{1.1892},{-0.037633},{0.32729}}
	Make/O/D/N=(1,3) matS={{0.17464},{-0.18671},{0.72579}}
End

Function MakeVectors()
	Make/O/D/N=(1,3) matP={{-0.43256},{-1.6656},{0.12533}}
	Make/O/D/N=(1,3) matQ={{0.28768},{-1.1465},{1.1909}}
	Make/O/D/N=(1,3) matR={{1.1892},{-0.037633},{0.32729}}
	Make/O/D/N=(1,3) matS={{0.17464},{-0.18671},{0.72579}}
End

Function DoCalcs()
	WAVE matP,matQ,matR,matS
	MatrixOp/O tempMat = matP - matQ
	MatrixSVD tempMat
	Make/O/D/N=(3,2) matN
	Wave M_VT
	matN = M_VT&#x5B;p]&#x5B;q+1]
	MatrixOp/O tempMat2 = (matR - matP)
	MatrixMultiply tempMat2, matN
	Wave M_product
	Duplicate/O M_product, mat_r
	MatrixOp/O tempMat2 = (matS - matP)
	MatrixMultiply tempMat2, matN
	Duplicate/O M_product, mat_s
	Make/O/D/N=(2,2) MatUnit
	matUnit = {{0,1},{-1,0}}
	MatrixOp/O tempMat2 = (mat_s - mat_r)
	MatrixMultiply tempMat2,MatUnit
	Duplicate/O M_Product, Mat_n
	Variable nn
	nn = norm(mat_n)
	MatrixOP/O new_n = mat_n / nn
	//new_n is now a vector with unit length
	Variable dd
	dd = MatrixDot(new_n,mat_r)
	//print dd
	//dd = MatrixDot(new_n,mat_s)
	//print dd
	dd = abs(dd)
	// now find v
	// v = dot(s-r,d*n-r)/dot(s-r,s-r)
	variable vv
	MatrixOp/O mat_s_r = mat_s - mat_r
	MatrixOp/O tempmat2 = dd * mat_n - mat_r
	vv = MatrixDot(mat_s_r,tempmat2) / MatrixDot(mat_s_r,mat_s_r)
	//print vv
	//because vv &amp;amp;gt; 1 then closest post is s (because rs(1) = s) and therefore closest point on RS to infinite line PQ is S
	//what about the point on PQ is this also outside the segment?
	// u = (Q-P)'\((S - (S*N)*N') - P)'
	variable uu
	MatrixOp/O matQ_P = matQ - matP
	MatrixTranspose matQ_P
	//MatrixOP/O tempMat2 = ((matS - (matS * matN) * MatrixTranspose(MatN)) - MatrixTranspose(matP))
	Duplicate/O MatN, matNprime
	MatrixTranspose matNprime
	MatrixMultiply matS, matN
	Duplicate/O M_Product, matSN
	MatrixMultiply M_Product, matNprime
	MatrixOP/O tempMat2 = ((matS - M_product) - matP)
	MatrixTranspose tempMat2
	MatrixLLS matQ_P tempMat2
	Wave M_B
	uu = M_B&#x5B;0]
	// find point on PQ that is closest to RS
	// P + u*(Q-P)
	MatrixOp/O matQ_P = matQ - matP
	MatrixOp/O matPoint = MatP + (uu * MatQ_P)
	MatrixOP/O distpoint = matPoint - matS
	Variable dist
	dist = norm(distpoint)
	Print dist
End
</pre>
The sticking points were finding the Igor equivalents of
<ul>
<li>null()</li>
<li>norm()</li>
<li>dot()</li>
<li>\ otherwise known as mldivide</li>
</ul>
Which are:
<ul>
<li>MatrixSVD (answer is in the final two columns of wave M_VT)</li>
<li>norm()</li>
<li>MatrixDot()</li>
<li>MatrixLLS</li>
</ul>
MatrixLLS wouldn&#8217;t accept a mix of single-precision and double-precision waves, so this needed to be factored into the code.
As you can see, the Igor code is much longer. Overall, I think MATLAB handles Matrix Math better than Igor. It is certainly easier to write. I suspect that there are a series of Igor operations that can do what I am trying to do here, but this was an exercise in direct porting.
More work is needed to condense this down and also deal with every possible case. Then it needs to be incorporated into the bigger program. SIMPLE! Anyway, hope this helps somebody.
&#8212;
The post title is taken from the band Adventures In Stereo.

An occasional series in esoteric programming issues.

As part of a larger analysis project I needed to implement a short
program to determine the closest distance of two line segments in 3D
space. This will be used to sort out which segments to compare... like I
say, part of a bigger project. The best method to do this is to find the
closest distance one segment is to the other when the other one is
represented as an infinite line. You can then check if that point is
beyond the segment if it is you use the limits of the segment to
calculate the distance. There's a discussion on stackoverflow
[here](http://stackoverflow.com/questions/627563/calculating-the-shortest-distance-between-two-lines-line-segments-in-3d).
The solutions point to [one in
C++](http://geomalgorithms.com/a07-_distance.html) and one in
[MATLAB](http://stackoverflow.com/a/702174). The C++ version is easiest
to port to Igor due to the similarity of languages, but the explanation
of the MATLAB code was more approachable. So I ported that to Igor to
figure out how it works.

The MATLAB version is:

``` {.brush: .matlabkey; .title: .; .notranslate title=""}
>> P = [-0.43256      -1.6656      0.12533]
P =
   -0.4326   -1.6656    0.1253
>> Q = [0.28768      -1.1465       1.1909]
Q =
    0.2877   -1.1465    1.1909
>> R = [1.1892    -0.037633      0.32729]
R =
    1.1892   -0.0376    0.3273
>> S = [0.17464     -0.18671      0.72579]
S =
    0.1746   -0.1867    0.7258
>> N = null(P-Q)
N =
   -0.3743   -0.7683
    0.9078   -0.1893
   -0.1893    0.6115
>> r = (R-P)*N
r =
    0.8327   -1.4306
>> s = (S-P)*N
s =
    1.0016   -0.3792
>> n = (s - r)*[0 -1;1 0];
>> n = n/norm(n);
>> n
n =
    0.9873   -0.1587
>> d = dot(n,r)
d =
    1.0491
>> d = dot(n,s)
d =
    1.0491
>> v = dot(s-r,d*n-r)/dot(s-r,s-r)
v =
    1.2024
>> u = (Q-P)'\((S - (S*N)*N') - P)'
u =
    0.9590
>> P + u*(Q-P)
ans =
    0.2582   -1.1678    1.1473
>> norm(P + u*(Q-P) - S)
ans =
    1.0710
```

and in IgorPro:

``` {.brush: .plain; .title: .; .notranslate title=""}
Function MakeVectors()
    Make/O/D/N=(1,3) matP={{-0.43256},{-1.6656},{0.12533}}
    Make/O/D/N=(1,3) matQ={{0.28768},{-1.1465},{1.1909}}
    Make/O/D/N=(1,3) matR={{1.1892},{-0.037633},{0.32729}}
    Make/O/D/N=(1,3) matS={{0.17464},{-0.18671},{0.72579}}
End

Function MakeVectors()
    Make/O/D/N=(1,3) matP={{-0.43256},{-1.6656},{0.12533}}
    Make/O/D/N=(1,3) matQ={{0.28768},{-1.1465},{1.1909}}
    Make/O/D/N=(1,3) matR={{1.1892},{-0.037633},{0.32729}}
    Make/O/D/N=(1,3) matS={{0.17464},{-0.18671},{0.72579}}
End

Function DoCalcs()
    WAVE matP,matQ,matR,matS
    MatrixOp/O tempMat = matP - matQ
    MatrixSVD tempMat
    Make/O/D/N=(3,2) matN
    Wave M_VT
    matN = M_VT[p][q+1]
    MatrixOp/O tempMat2 = (matR - matP)
    MatrixMultiply tempMat2, matN
    Wave M_product
    Duplicate/O M_product, mat_r
    MatrixOp/O tempMat2 = (matS - matP)
    MatrixMultiply tempMat2, matN
    Duplicate/O M_product, mat_s
    Make/O/D/N=(2,2) MatUnit
    matUnit = {{0,1},{-1,0}}
    MatrixOp/O tempMat2 = (mat_s - mat_r)
    MatrixMultiply tempMat2,MatUnit
    Duplicate/O M_Product, Mat_n
    Variable nn
    nn = norm(mat_n)
    MatrixOP/O new_n = mat_n / nn
    //new_n is now a vector with unit length
    Variable dd
    dd = MatrixDot(new_n,mat_r)
    //print dd
    //dd = MatrixDot(new_n,mat_s)
    //print dd
    dd = abs(dd)
    // now find v
    // v = dot(s-r,d*n-r)/dot(s-r,s-r)
    variable vv
    MatrixOp/O mat_s_r = mat_s - mat_r
    MatrixOp/O tempmat2 = dd * mat_n - mat_r
    vv = MatrixDot(mat_s_r,tempmat2) / MatrixDot(mat_s_r,mat_s_r)
    //print vv
    //because vv &amp;gt; 1 then closest post is s (because rs(1) = s) and therefore closest point on RS to infinite line PQ is S
    //what about the point on PQ is this also outside the segment?
    // u = (Q-P)'\((S - (S*N)*N') - P)'
    variable uu
    MatrixOp/O matQ_P = matQ - matP
    MatrixTranspose matQ_P
    //MatrixOP/O tempMat2 = ((matS - (matS * matN) * MatrixTranspose(MatN)) - MatrixTranspose(matP))
    Duplicate/O MatN, matNprime
    MatrixTranspose matNprime
    MatrixMultiply matS, matN
    Duplicate/O M_Product, matSN
    MatrixMultiply M_Product, matNprime
    MatrixOP/O tempMat2 = ((matS - M_product) - matP)
    MatrixTranspose tempMat2
    MatrixLLS matQ_P tempMat2
    Wave M_B
    uu = M_B[0]
    // find point on PQ that is closest to RS
    // P + u*(Q-P)
    MatrixOp/O matQ_P = matQ - matP
    MatrixOp/O matPoint = MatP + (uu * MatQ_P)
    MatrixOP/O distpoint = matPoint - matS
    Variable dist
    dist = norm(distpoint)
    Print dist
End
```

The sticking points were finding the Igor equivalents of

- null()
- norm()
- dot()
- \\ otherwise known as mldivide

Which are:

- MatrixSVD (answer is in the final two columns of wave M_VT)
- norm()
- MatrixDot()
- MatrixLLS

MatrixLLS wouldn't accept a mix of single-precision and double-precision
waves, so this needed to be factored into the code.

As you can see, the Igor code is much longer. Overall, I think MATLAB
handles Matrix Math better than Igor. It is certainly easier to write. I
suspect that there are a series of Igor operations that can do what I am
trying to do here, but this was an exercise in direct porting.

More work is needed to condense this down and also deal with every
possible case. Then it needs to be incorporated into the bigger program.
SIMPLE! Anyway, hope this helps somebody.

---

The post title is taken from the band Adventures In Stereo.


An occasional series in esoteric programming issues. As part of a larger analysis project I needed to implement a short program to determine the closest distance of two line segments in 3D space. This will be used to sort out which segments to compare… like I say, part of a bigger project. The best method to do this is to find the closest distance one segment is to the other when the other one is represented as an infinite line.


Rogue Scholar Beiträge

Adventures in code