Note that there are some explanatory texts on larger screens.

plurals
  1. PO
    text
    copied!<p>You can use scipy.spatial.cKDTree to speedup this kind of calculation.</p> <p>Here is the code:</p> <pre><code>import time import numpy as np #### create some sample data #### np.random.seed(1) V_NUM = 6000 P_NUM = 8000 size = 0.1 vertices = np.random.rand(V_NUM, 3) points = np.random.rand(P_NUM, 3) numpoints = np.zeros(V_NUM, np.int32) #### brute force #### start = time.clock() for i in xrange(len(vertices)): cube=((vertices[i,0]&lt;= points[:,0]) &amp; (points[:,0]&lt;(vertices[i,0]+size)) &amp; (vertices[i,1]&lt;= points[:,1]) &amp; (points[:,1] &lt; (vertices[i,1]+size)) &amp; (vertices[i,2]&lt;= points[:,2]) &amp; (points[:,2] &lt; (vertices[i,2]+size)) ) numpoints[i]=len(points[cube]) print time.clock() - start #### KDTree #### from scipy.spatial import cKDTree center_vertices = vertices + [[size/2, size/2, size/2]] start = time.clock() tree_points = cKDTree(points) _, result = tree_points.query(center_vertices, k=100, p = np.inf, distance_upper_bound=size/2) numpoints2 = np.zeros(V_NUM, np.int32) for i, neighbors in enumerate(result): numpoints2[i] = np.sum(neighbors!=P_NUM) print time.clock() - start print np.all(numpoints == numpoints2) </code></pre> <ul> <li>Change the cube corner position to center position first. </li> </ul> <p><code>center_vertices = vertices + [[size/2, size/2, size/2]]</code></p> <ul> <li>Create a cKDTree from points</li> </ul> <p><code>tree_points = cKDTree(points)</code></p> <ul> <li>Do query, k is the the number of nearest neighbors to return, p=np.inf means maximum-coordinate-difference distance, distance_upper_bound is the max distance.</li> </ul> <p><code>_, result = tree_points.query(center_vertices, k=100, p = np.inf, distance_upper_bound=size/2)</code></p> <p>The output is:</p> <pre><code>2.04113164434 0.11087783696 True </code></pre> <p>If there are more then 100 points in a cube, you can check this by <code>neighbors[-1] == P_NUM</code> in the for loop, and do a k=1000 query for these vertices.</p>
 

Querying!

 
Guidance

SQuiL has stopped working due to an internal error.

If you are curious you may find further information in the browser console, which is accessible through the devtools (F12).

Reload