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. The solutions point to one in C++ and one in MATLAB. 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:

>> 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:

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 &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.