Find all possible distances between two sets of xy coords

Share icon

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.

February 28, 2018 at 12:31 am - Permalink

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.

February 28, 2018 at 04:35 am - Permalink

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.

February 28, 2018 at 08:33 am - Permalink

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.

February 28, 2018 at 09:08 am - Permalink

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.

February 28, 2018 at 10:42 am - Permalink

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

February 28, 2018 at 11:36 am - Permalink

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.

February 28, 2018 at 12:24 pm - Permalink