Given a cloud of points, what would be the smartest way to determine the two points that are most widely separated?
It is possible to loop through all points and calculate the distances between points so as to get the maximum, it would require n x n iterations. I wish to know if there exists a class that can do that, preferably not Python exclusive.
This is a classic computational geometry problem with some smart solutions. I would expect that you do not have to implement one, but you can find free, open- source implementations in Python and C++.
vtkSphere offers to compute a bounding sphere. The problem is that it requires a hints[2] field, which is a guess of the 2 most distant points. Random pairs give different results, so getting the exact pair of points is my purpose. Iāll delve further into this slowly.
Implementation in VTK seems to be a standard one (same as in three.js) that computes a bounding sphere with a simple method, which is not guaranteed to be minimal, but it is close. Hints are not required. If you donāt provide hints then the bounding box extreme points are used to initialize the computation.
If you only need an approximate binding sphere (that can be a few percent larger than needed) then this should work well. If you need the optimal sphere then it will require a more complex, more costly computation.
What are you planning to use this bounding sphere for?
Itās intended to add 1 more option to this module.
It needs not have an exact precision. However, the āhints[2]ā parameter can be passed as ānullptrā in C++, while passing āNoneā fails in Python. Passing arbitrary values like ā[0, 0]ā is meaningless. How to handle this situation in Python?
Unfortunately, I needed to pass a āNoneā parameter that would mean ānullptrā, but the python wrapper always wants to pass āsomethingā else than 0x0 to the underlying C++ function.
It would have been nice if this C++ function signature defined āhints[2] = nullptrā since it is optional. I suppose it should have worked.
Anyway, Iām baking a solution that should not be resource intensive.
Agreed. You can send a merge request to VTK with this simple change. It should go through review really quickly. You can also @mention David Gobbi in the merge request, asking him about the VTK null pointer syntax in Python (maybe there is a way to specify a null pointer just not documented well).
I came up to this function to determine the two most distant points in a cloud of points. Itās decently fast, 0.06 s for 10ā¶ points, 1.01 s for 10ā· points, 11.69 s for 10āø points and bumps to 135.67 s for 10ā¹ points.
The logic seems OK to me, dropping it here for any purpose and remaining open for comments.
def findMostDistantPoints(points : vtk.vtkPoints):
pointData = points.GetData()
numberOfPoints = points.GetNumberOfPoints()
# Sort the points on each axis, get the bounding points.
sorter = vtk.vtkSortDataArray()
sorter.SortArrayByComponent(pointData, 0, 0)
xMin = pointData.GetTuple3(0)
xMax = pointData.GetTuple3(numberOfPoints - 1)
sorter.SortArrayByComponent(pointData, 1, 0)
yMin = pointData.GetTuple3(0)
yMax = pointData.GetTuple3(numberOfPoints - 1)
sorter.SortArrayByComponent(pointData, 2, 0)
zMin = pointData.GetTuple3(0)
zMax = pointData.GetTuple3(numberOfPoints - 1)
boundingPoints = vtk.vtkPoints()
boundingPoints.InsertNextPoint(xMin)
boundingPoints.InsertNextPoint(xMax)
boundingPoints.InsertNextPoint(yMin)
boundingPoints.InsertNextPoint(yMax)
boundingPoints.InsertNextPoint(zMin)
boundingPoints.InsertNextPoint(zMax)
# Calculate the distances between each pair of points: 6 x 6 iterations.
distances = vtk.vtkDoubleArray()
distances.SetNumberOfComponents(7)
numberOfBoundingPoints = boundingPoints.GetNumberOfPoints()
for i in range(numberOfBoundingPoints):
referencePoint = boundingPoints.GetPoint(i)
for j in range(numberOfBoundingPoints):
nextPoint = boundingPoints.GetPoint(j)
distance2 = vtk.vtkMath.Distance2BetweenPoints(referencePoint, nextPoint)
distances.InsertNextTuple([referencePoint[0], referencePoint[1], referencePoint[2], nextPoint[0], nextPoint[1], nextPoint[2], distance2])
# Sort the bounding points in descending order by distance squared.
distanceSorter = vtk.vtkSortDataArray()
distanceSorter.SortArrayByComponent(distances, 6, 1)
# Get the most distant points.
mostDistantTuple = distances.GetTuple(0)
p1 = [mostDistantTuple[0], mostDistantTuple[1], mostDistantTuple[2]]
p2 = [mostDistantTuple[3], mostDistantTuple[4], mostDistantTuple[5]]
# Return as vtkPoints.
result = vtk.vtkPoints()
result.InsertNextPoint(p1)
result.InsertNextPoint(p2)
return result
Computation is quite fast and the method is simple, but I donāt think it gives you the two farthest points (when the two farthest points are near the corners of the bounding box, in two opposite corners).
Instead of reinventing the wheel, you can use a free, open-source, restriction-free implementation like this.
How does it work? Have you you implemented a known algorithm, or you invented something new? If it is something new, how it compares to known algorithms?
On test data, it seems to be accurate. Iāll compare when VTK 8fcf3271ab is available in Slicerās VTK. I just wanted to propose another possible solution to a question I opened.