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 > 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.
The C++ version (mentioned above) was much easier to rewrite… only taking about 30 min to do. This version was more readily usable.