# Find all possible distances between two sets of xy coords

I have two 2D waves of xy coords with a different number of points and need to find all possible distances between them. I feel like there should be a MatrixOp command for this that I am missing. I have thousands of pairs to process so speed is important.

Currently I do this:

make/o/n=(10,2) aa=p*(q+2) // make some integer coords for illustration
make/o/n=(20,2) bb
•bb[][0] += 3
•bb[][1] += 4 // second set is 5 units away
make/o/n=(dimsize(aa,0),dimsize(bb,0),2) cc // to contain the difference between each point, x in layer 0, y in layer 1
matrixtranspose bb
•cc[][][0] = aa[p][0] - bb[0][p]
•cc[][][1] = aa[p][1] - bb[1][p] // subtract each point in bb from each point in aa
matrixop/o dd = cc * cc
matrixop/o ee = sumbeams(dd)
matrixop/o ff = sqrt(ee) // 2D wave of all distances

I can't do the last three MatrixOps as a compound expression because of sumbeams.

Is there a simple command to do this? Or any other way to simplify these steps? Any input would be appreciated.
Elsewhere in my code I turn these coordinate sets into image masks using `ImageBoundaryToMask` and I looked at image operations to find the distances but couldn't figure out a good way to do this.
My first thought is that you are working on a clustering problem. If that is the case, you may want to look at FPClustering or KMeans.

Otherwise, here is my way of using matrixop. I broke it down to individual instructions and it is past midnight here so you may want to check it twice.
Function test()
// make sample data
make/O/n=10 ax=p,ay=p
make/o/n=8 bx=10*p,by=10*p

// create matrices that are 10x8:
Matrixop/o/free matAx=colRepeat(ax,8)
Matrixop/o/free matAy=colRepeat(ay,8)
Matrixop/o/free matBx=rowRepeat(bx,10)
Matrixop/o/free matBy=rowRepeat(by,10)
Matrixop/o/free distanceX=matAx-matBx
Matrixop/o/free distanceY=matAy-matBy
MatrixOP/o distances=sqrt(distanceX*distanceX+distanceY*distanceY)
End

Now you have to pick the relevant data (array) out of the full matrix. It is not the most efficient but it is probably worth trying.

I hope this helps,

A.G.
WaveMetrics, Inc.
Thank you for the late night reply. I like your approach. It's more readable than mine. I'm glad I'm not missing an obvious command.

I'm not looking at clusters right now. I actually just want to find the closest distance between the two sets of points.
sjr51 wrote:
Thank you for the late night reply. I like your approach. It's more readable than mine. I'm glad I'm not missing an obvious command.

I'm not looking at clusters right now. I actually just want to find the closest distance between the two sets of points.

Look at FPClustering. I think that if you set it correctly you could get the closest distance.
A.G.'s solution is clever and elegant. His example however is special in that it gives an uninteresting solution with distance 0. Here is a slightly more general case, with a simple means to return the minimum distance:
Function test()
// make sample data
make/O/n=10 ax=p+2.5,ay=p+3.75
make/o/n=8 bx=10*p,by=10*p
// create matrices that are 10x8:
Matrixop/o/free matAx=colRepeat(ax,8)
Matrixop/o/free matAy=colRepeat(ay,8)
Matrixop/o/free matBx=rowRepeat(bx,10)
Matrixop/o/free matBy=rowRepeat(by,10)
Matrixop/o/free distanceX=matAx-matBx
Matrixop/o/free distanceY=matAy-matBy
MatrixOP/o distances=sqrt(distanceX*distanceX+distanceY*distanceY)

ImageStats distances
return V_min
End

print test()
0.901388

ImageStats can also provide the matrix location of the minimum distance at V_minColLoc and V_minRowLoc. You will have to figure out whether or how to search for other equal minimum distances in your data set.
Just as an update:

FPClustering can't work here because it normalizes each column of the position data. It is also difficult to make it output distances that are not clustered. I will look into adding such an option for IP8.

A.G.
WaveMetrics, Inc.
If you have access to the beta, after the nightly build, try the following:
Function testInIP8()
// Make sample data using two-column waves.  For higher dimensions use cols=nDims.
make/O/n=(10,2) aa
aa[][0]=p+2.5
aa[][1]=p+3.75
make/o/n=(8,2) bb
bb[][0]=10*p
bb[][1]=10*p

// Combine the two distance groups into a single positions wave:
Concatenate/NP=0 {aa,bb}, positions
FPClustering/DSO positions                          // Requires IP8
Wave M_DistanceMap
// the output M_DistanceMap is an upper triangular matrix containing the distances.
// since you do not care about the distances between positions internal to {aa} or {bb}
// extract the relevant matrix:
Variable colStart=DimSize(aa,0)
Variable colEnd=colStart+DimSize(bb,0)-1
Variable rowsStart=0
Variable rowsEnd=colStart-1
// you can extract the relevant part by:
MatrixOP/o/FREE myDistance=subRange(M_DistanceMap,rowsStart,rowsEnd,colStart,colEnd)
return WaveMin(myDistance)
End
I do have access and will try. Thank you for implementing it.

@s.r.chinn thank you for the input. I am using exactly this in my code right now.