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