Is there a fast way to get the minimum distance between two vertices in a geometry?

I need a method to get the minimum distance between all the vertices in a geometry and write a one with two "for" loops. It works well when the geometry is simple, whereas it takes too much time if there are thousands of vertices.

Thus I want some help in this: if there is a simpler way provided or not? Or how can I do it with (O) complexity. Thanks a lot for providing information. Looking forward to replies.

jelle's picture

what if you store the vertices of the 2nd shape as a kd-tree? that would speed up things quite a bit...

Marco Nawijn's picture

I can confirm that a kdtree structure is very well suited for this case. In python you can find an implementation in both BioPython and scipy. I have generated a small wrapper class to make working with the kdtrees a little easier. You can find the code below. Note that it works easiest when you also have a Point class representing the point in 3D space.

class KDTree(object):
'''Represents a generic KDTree data structure. The current class is a small
wrapper around the Biopython KDTree class. It makes the interface a little
more convenient when using the KDTree class with higher level objects.

'''

def __init__(self, dim=3, bucket_size=1):
'''Initializes a new tree and supporting data structures.

'''
self.dim = dim

self._tree = _kdtree.KDTree(dim, bucket_size)

self._domain = {}
self.radii = None

self._idx_pid_map = {}

self.closest_points = None
self.closest_radii = None

@property
def domain(self): return self._domain.values()

@domain.setter
def domain(self, domain):
'''Property setter for domain.'''

context = 'KDTree._setDomain()'

try:
M = len(domain) # Number of rows
except TypeError:
error('Domain should be a list or a tuple.', context)

if M==0:
error('No points in domain.', context)

N = self.dim # Number of columns

crds = empty((M,N),dtype='float32')

i = 0

for p in domain:
crds[i,:] = p.crds

self._idx_pid_map[i] = p.pid

self._domain[p.pid] = p

i+=1

self._tree.set_coords(crds)

def search(self, target, radius):
'''Searches the domain for points that are located within the distance
specified by radius to the target point or coordinate.

The target can be either a Point instance or a cartesian coordinate.

All points and coordinates are treated as 3 dimensional. This implies
that additional zero entries are automatically inserted in the coordinate
arrays.

Typical usage:

>> t = KDTree()
>> t.domain = domain # domain is a list of Point instances
>> target = Point(1, (5.0, 5.0, 0.0))
>> t.search(target, 0.5)
Point 383, dim = 3, 4.778950, 4.950348, 0.000000, ndata = 0
Point 423, dim = 3, 5.075207, 5.242366, 0.000000, ndata = 0
Point 600, dim = 3, 5.160043, 4.737402, 0.000000, ndata = 0

'''

# Make sure we start with an empty result data sets
self.radii = {}
self.closest_points = []
self.closest_radii = []

if isinstance(target, Point):
target = target.crds

else:
if len(target)==1:
if isinstance(target, tuple):
target = list(target)
target.extend((0.0, 0.0))
if len(target)==2:
if isinstance(target, tuple):
target = list(target)
target.append(0.0)
if len(target)>3:
error('Invalid number of components in target coordinates')

_c = array(target, dtype='float32')
self._tree.search(_c, radius)

idxs = self._tree.get_indices()
radii = self._tree.get_radii()

result = []
for idx, r in izip(idxs, radii):
pid = self._idx_pid_map[idx]
result.append(self._domain[pid])
self.radii[pid] = r

return result

def closest(self, target, radius, npoints=1):
'''Obtain the "npoints" closest points to the target point. By default,
only a single point will be returned. If more points are required a list
of these points will be returned.

The number of points N can be one of the following:
N = 0
if there are no points within the specified radius

N = 1
if only one point is requested or if there is only 1 point within
the search radius

N = more than 1 but smaller than npoints
if there are less points within the radius than the number of
requested points

N>=npoints
if there are more points within the radius than the number of
requested points

Typical usage:

>> t = KDTree()
>> t.domain = domain # domain is a list of Point instances
>> target = Point(1, (5.0, 5.0, 0.0))
>> t.closest(target, 0.5, npoints=3)
Point 383, dim = 3, 4.778950, 4.950348, 0.000000, ndata = 0
Point 423, dim = 3, 5.075207, 5.242366, 0.000000, ndata = 0
Point 600, dim = 3, 5.160043, 4.737402, 0.000000, ndata = 0
'''

self.search(target, radius)

result = []
for pid, r in self.radii.items():
result.append((r,pid))

nresult = len(result)
if nresult>0:
result.sort()
idx=0
for r in result:
if idx==npoints: break
radius, pid = r
point = self._domain[pid]
self.closest_radii.append(radius)
self.closest_points.append(point)
idx+=1

if len(self.closest_points)==1 or npoints==1:
return self.closest_points[0]
else:
return self.closest_points
else: return None

Marco Nawijn's picture

Hmm...indentation is screwed while posting. If you want the code drop me a PM.

Regards,

Marco

jelle's picture

cool! thanks for sharing that wrapper Marco...
perhaps you can add the .py file as an attachment?

cheers,
-jelle

Marco Nawijn's picture

Here it is. I included the definition of the Point class for your convenience ;). It should be trivial to change the code so it works with OCC vertices. I can take a look if someone is interested and add it to pythonocc. One thing I am currently looking at is changing the Biopython kdtree implementation for the one provide by scipy.

Regards,

Marco

Attachments: