Speeds up the selection of sets of three nodes that form triangles with a given minimum and maximum side length
I have a list of 60 points, from which I create all possible combinations of size 3. I am trying to determine if any of the three points are within certain distance parameters, that is, not too close to each other and not too far apart (say no point can be less than 10 units to the nearest nearest point, and no point can be more than 100 units from the farthest point).
I have some code to do it here:
def setPalmCombinations(self):
# Testing confirms that the range of distances is about 5 - 10 for
# minimum distance between points, and anywhere from 150 - 200 for max
# distance between points. 10 and 100 should be reasonable cutoffs.
minDistance = 10
maxDistance = 100
startTime = time.time()
iter = itertools.combinations(self.handContour, 3)
def distanceCheck(points):
A = points[0][0]
B = points[1][0]
C = points[2][0]
AB = np.linalg.norm(A - B)
if not(minDistance < AB < maxDistance): return None
BC = np.linalg.norm(B - C)
if not(minDistance < BC < maxDistance): return None
CA = np.linalg.norm(C - A)
if not(minDistance < CA < maxDistance): return None
return np.array([A, B, C])
a = [distanceCheck(i) for i in iter]
print time.time() - startTime
However, just creating this list of possible combinations takes about 0.4 ... 0.5 seconds, which is too slow. Is there a way I can optimize this calculation (or fundamentally rewrite the algorithm so that it is not just a brute force search) so that I can shorten the time to compute in real time (i.e. 30 times per second).
I was thinking about creating a list of distances between each point, sorting it, and finding the points with distances to the center of the clipping range, but I think this might actually be slower than what I am doing now.
I'm also trying to figure out some smarter algorithm, given the fact that the points are initially stored in the order of the contour definition - I should be able to eliminate some points before and after each point in the list when making combinations ... however when I tried to generate combinations using this fact with Python's default for loops, it turned out to be slower than just using itertools and generated everything.
Any advice is appreciated.
Thank!
source to share
I would first calculate the pairwise distances between all points, then filter the ones that fit your specifications ( dist_min
and dist_max
), and then check to see if the combination of the three criteria of the nodes corresponds to (because I just checked in pairs - once for speed).
import itertools
import numpy as np
from scipy.spatial.distance import pdist, squareform
points = np.random.rand(100,2) # could be nodes in 3 dimensions too, but you specify a handpalm, so I'm guessing your contour is merely 2D
def palm_combinations(points, dist_min=.0, dist_max=.05):
dists = pdist(points)
# This matrix lists the distance between each pair of nodes
# (represented by their column and row index, so that the mat[i,j]
# represents the distance between points[i,:], points[j,:])
mat = squareform(dists)
mask = (mat > dist_min) & (mat < dist_max)
two_good_legs = {} # store all combinations of point indices where "two legs" of the triangles have a distance that matches the spec.
myinds = []
for rowind, row in enumerate(mask):
relevant = row[rowind+1:] # distance to a point itself is always zero, hence exclude from list
inds = np.where(relevant)[0] + (rowind+1)
myinds.append(inds)
two_good_legs[rowind] = list(itertools.combinations(inds, 2))
triangles = []
# The only thing left to do now is to check if the third leg norm is also within spec.
for k,v in two_good_legs.items():
for tup in v:
if tup[1] in myinds[tup[0]]:
triangles.append((k, tup[0], tup[1]))
return triangles
triangles
will display all indices at points that match your requirements. If you want to go back to the actual coordinates, just get points[triangles,:]
. You can visualize all of this with:
import maptlotlib.pyplot as plt
plt.triplot(points[:,0], points[:,1], triangles)
plt.plot(points[:,0], points[:,1], 'ko ')
Temporal differences between your implementation (modified to exclude unknown [0]
indexing, see my comment) show this:
In [12]: %timeit my_implementation(points, .1, .15)
1 loops, best of 3: 1.58 ms per loop
In [13]: %timeit your_implementation(points, .1, .15)
1 loops, best of 3: 2.32 s per loop
With points shaped (100,2)
so that the speed increased by 1470x on my machine.
You can increase the speed a little if you get rid of the call squareform
, which is a bottleneck according to the profiler and is not absolutely necessary, but makes the code more readable ("Read counts." Cfr. Python Zen, Tim Peters). To do this, you will make the following changes:
def palm_combinations(points, dist_min=.0, dist_max=.05):
dists = pdist(points)
mask = (dists > dist_min) & (dists < dist_max)
two_good_legs, myinds, triangles = {}, [], []
s = points.shape[0]
ind_u1, ind_un = 0, s - 1
for rowind in xrange(s-1):
relevant = mask[ind_u1:ind_un]
ind_u1 = ind_un
ind_un += s - 1*(rowind+2)
inds = np.where(relevant)[0] + (rowind+1)
myinds.append(inds)
two_good_legs[rowind] = list(itertools.combinations(inds, 2))
for k,v in two_good_legs.items():
for tup in v:
if tup[1] in myinds[tup[0]]:
triangles.append((k, tup[0], tup[1]))
return triangles
However, the increase in speed will not be noticeable for small datasets.
source to share
You can speed up using numpy. Assuming comb
- these are your possible combinations already calculated, you can easily vectorize the distance calculations:
>>> import numpy as np
>>> comb = np.array(list(comb))
>>> dists = comb - np.roll(comb, -1)
# L2 norm
>>> dists = dists**2
dists will array Nx3
containing AB
BC
and CA
distance, respectively.
Customize the mask with your conditions:
>>> mask = (dists > minDist) & (dists < maxDist)
Find lines where all distances satisfy the condition:
>>> rows = np.all(mask, axis=1)
which satisfy the conditions:
>>> comb[rows]
UPDATE: just noticed that this is quite slow for large N. Here is a function that contains all of the above steps:
>>> def palm_combinations(points, min_dist=10, max_dist=100):
combs = np.array(list(itertools.combinations(points, 3)))
dists = (combs - np.roll(combs, -1))**2
mask = (dists > min_dist) & (dists < max_dist)
return combs[np.all(mask, axis=1)]
For sample data of 100 items:
>>> a = np.random.rand(100) * 500
>>> %timeit palm_combinations(a)
10 loops, best of 3: 79.5 ms per loop
The bottleneck of the function is the generation of all combinations:
>>> %timeit combs = np.array(list(itertools.combinations(a, 3)))
10 loops, best of 3: 69.2 ms per loop
UPDATE 2: It can be sped up with a more sophisticated approach. Using np.fromiter
:
# Define combinations data type (3 float each combination)
>>> dt = np.dtype('f,f,f')
# Fast numpy array conversion
>>> combs = np.fromiter(itertools.combinations(points, 3), dt)
The function will look something like this:
>>> def palm_combinations(points, min_dist=10, max_dist=100):
dt = np.dtype('f,f,f')
combs = np.fromiter(itertools.combinations(points, 3), dt)
combs = combs.view('f').reshape(-1,3)
dists = (combs - np.roll(combs, -1))**2
mask = (dists > min_dist) & (dists < max_dist)
return combs[np.all(mask, axis=1)]
And the speed test:
>>> %timeit palm_combinations(a)
10 loops, best of 3: 29.8 ms per loop
For N = 100
source to share
If the number of contour points is very large, you can try the heuristic approach. This will take three random points, (check the minimum and maximum distances), calculate the radius of the corresponding circle, compare with the previously calculated radius, repeat. You would have to figure in some kind of stopping criteria, such as the maximum number of iterations or no better radius found in the last x iterations, or the growth rate of the radii below x percent, or some combination of these.
Also, in regards to your clarifying comment on your question: Either your terminology or your approach is off. There is no guarantee that the circle passing through a contour point is inside the convex hull of the contour points. Imagine, for example, an extended arm, the shape will have concave segments, if all three points are within one such segment, the circle will lie outside the arm. I think you mean a circle inscribed in the convex hull of every three point set.
source to share