Sun, 08/21/2011 - 05:55

Forums:

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.

Sun, 08/21/2011 - 13:58

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

Mon, 08/22/2011 - 16:27

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

Mon, 08/22/2011 - 16:28

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

Regards,

Marco

Mon, 08/22/2011 - 16:57

cool! thanks for sharing that wrapper Marco...

perhaps you can add the .py file as an attachment?

cheers,

-jelle

Mon, 08/22/2011 - 17:10

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