pep8 compliance
authorMichel Anders <michel.anders@inter.nl.net>
Fri, 26 Jul 2013 19:24:25 +0000 (19:24 +0000)
committerMichel Anders <michel.anders@inter.nl.net>
Fri, 26 Jul 2013 19:24:25 +0000 (19:24 +0000)
add_mesh_space_tree/__init__.py
add_mesh_space_tree/kdtree.py
add_mesh_space_tree/sca.py
add_mesh_space_tree/simplefork.py
add_mesh_space_tree/timer.py

index cd565f4..0637c41 100644 (file)
@@ -24,7 +24,7 @@
 bl_info = {
     "name": "SCA Tree Generator",
     "author": "michel anders (varkenvarken)",
-    "version": (0, 1, 0),
+    "version": (0, 1, 1),
     "blender": (2, 66, 0),
     "location": "View3D > Add > Mesh",
     "description": "Adds a tree created with the space colonization algorithm starting at the 3D cursor",
index 517605f..648f714 100644 (file)
 #
 # ##### END GPL LICENSE BLOCK #####
 
+# <pep8 compliant>
+
 from copy import copy, deepcopy
 
+
 class Hyperrectangle:
-       '''an axis aligned bounding box of arbitrary dimension'''
-
-       def __init__(self, dim, min, max):
-               self.dim = dim
-               self.min = deepcopy(min) # min and max should never point to the same instance
-               self.max = deepcopy(max)
-
-       def extend(self,pos):
-               '''adapt the hyperectangle if necessary so it will contain pos.'''
-               for i in range(self.dim):
-                       if pos[i]<self.min[i]: self.min[i]=pos[i]
-                       elif pos[i]>self.max[i]: self.max[i]=pos[i]
-               
-       def distance_squared(self,pos):
-               '''return the distance squared to the nearest edge, or zero if pos lies within the hyperrectangle'''
-               result=0.0
-               for i in range(self.dim):
-                       if pos[i]<self.min[i]:
-                               result+=(pos[i ]-self.min[i ])**2
-                       elif pos[i]>self.max[i]:
-                               result+=(pos[i ]-self.max[i ])**2
-               return result
-
-       def __str__(self):
-               return "[(%d) %s:%s]"%(int(self.dim),str(self.min),str(self.max))
+    '''an axis aligned bounding box of arbitrary dimension'''
+
+    def __init__(self, dim, min, max):
+        self.dim = dim
+        self.min = deepcopy(min)  # min and max should never point to the same instance
+        self.max = deepcopy(max)
+
+    def extend(self, pos):
+        '''adapt the hyperectangle if necessary so it will contain pos.'''
+        for i in range(self.dim):
+            if pos[i] < self.min[i]:
+                self.min[i] = pos[i]
+            elif pos[i] > self.max[i]:
+                self.max[i] = pos[i]
+
+    def distance_squared(self, pos):
+        '''return the distance squared to the nearest edge, or zero if pos lies within the hyperrectangle'''
+        result = 0.0
+        for i in range(self.dim):
+            if pos[i] < self.min[i]:
+                result += (pos[i] - self.min[i]) ** 2
+            elif pos[i] > self.max[i]:
+                result += (pos[i] - self.max[i]) ** 2
+        return result
+
+    def __str__(self):
+        return "[(%d) %s:%s]" % (int(self.dim), str(self.min), str(self.max))
+
 
 class Node:
-       """implements a node in a kd-tree"""
-       
-       def __init__(self, pos, data=None):
-               self.pos = deepcopy(pos)
-               self.data = data
-               self.left = None
-               self.right = None
-               self.dim=len(pos)
-               self.dir = 0
-               self.count=0
-               self.level=0
-               self.rect=Hyperrectangle(self.dim,pos,pos)
-
-       def addleft(self, node):
-               self.left = node
-               self.rect.extend(node.pos)
-               node.level=self.level+1
-               node.dir=(self.dir+1)%self.dim
-               
-       def addright(self, node):
-               self.right = node
-               self.rect.extend(node.pos)
-               node.level=self.level+1
-               node.dir=(self.dir+1)%self.dim
-               
-       def distance_squared(self, pos):
-               d=self.pos-pos
-               return d.dot(d)
-
-       def _str(self,level):
-               s = '  '*level+str(self.dir)+' '+str(self.pos)+' '+str(self.rect)+'\n'
-               return s + ('' if self.left is None else 'L:'+self.left._str(level+1)) + ('' if self.right is None else 'R:'+self.right._str(level+1))
-
-       def __str__(self):
-               return self._str(0)
+    """implements a node in a kd-tree"""
+
+    def __init__(self, pos, data=None):
+        self.pos = deepcopy(pos)
+        self.data = data
+        self.left = None
+        self.right = None
+        self.dim = len(pos)
+        self.dir = 0
+        self.count = 0
+        self.level = 0
+        self.rect = Hyperrectangle(self.dim, pos, pos)
+
+    def addleft(self, node):
+        self.left = node
+        self.rect.extend(node.pos)
+        node.level = self.level + 1
+        node.dir = (self.dir + 1) % self.dim
+
+    def addright(self, node):
+        self.right = node
+        self.rect.extend(node.pos)
+        node.level = self.level + 1
+        node.dir = (self.dir + 1) % self.dim
+
+    def distance_squared(self, pos):
+        d = self.pos - pos
+        return d.dot(d)
+
+    def _str(self, level):
+        s = '  ' * level + str(self.dir) + ' ' + str(self.pos) + ' ' + str(self.rect) + '\n'
+        return s + ('' if self.left is None else 'L:' + self.left._str(level + 1)) + ('' if self.right is None else 'R:' + self.right._str(level + 1))
+
+    def __str__(self):
+        return self._str(0)
+
 
 class Tree:
-       """implements a kd-tree"""
-       
-       def __init__(self, dim):
-               self.root = None
-               self.nnearest=0 # number of nearest neighbor queries
-               self.count=0  # number of nodes visited
-               self.level=0 # deepest node level 
-       
-       def resetcounters(self):
-               self.nnearest=0 # number of nearest neighbor queries
-               self.count=0  # number of nodes visited
-               
-       def _insert(self, node, pos, data):
-               if pos[node.dir] < node.pos[node.dir]:
-                       if node.left is None:
-                               node.addleft(Node(pos, data))
-                               return node.left
-                       else:
-                               node.rect.extend(pos)
-                               return self._insert(node.left, pos, data)
-               else:
-                       if node.right is None:
-                               node.addright(Node(pos, data))
-                               return node.right
-                       else:
-                               node.rect.extend(pos)
-                               return self._insert(node.right, pos, data)
-
-       def insert(self, pos, data):
-               if self.root is None:
-                       self.root = Node(pos,data)
-                       self.level = self.root.level
-                       return self.root
-               else:
-                       node=self._insert(self.root, pos, data)
-                       if node.level > self.level : self.level = node.level
-                       return node
-
-       def _nearest(self, node, pos, checkempty, level=0):
-               
-               self.count+=1
-               
-               dir = node.dir
-               d = pos[dir] - node.pos[dir]
-               
-               result = node
-               distsq = None
-               if checkempty and (node.data is None):
-                       result = None
-               else:
-                       distsq = node.distance_squared(pos)
-
-               if d <= 0:
-                       neartree = node.left
-                       fartree = node.right
-               else:
-                       neartree = node.right
-                       fartree = node.left
-
-               if neartree is not None:
-                       nearnode, neardistsq = self._nearest(neartree,pos,checkempty,level+1)
-                       if (result is None) or (neardistsq is not None and neardistsq < distsq):
-                               result, distsq = nearnode, neardistsq
-               
-               if fartree is not None:
-                       if (result is None) or (fartree.rect.distance_squared(pos) < distsq):
-                               farnode, fardistsq = self._nearest(fartree,pos,checkempty,level+1)
-                               if (result is None) or (fardistsq is not None and fardistsq < distsq):
-                                       result, distsq = farnode, fardistsq
-                       
-               return result, distsq
-
-       def nearest(self, pos, checkempty=False):
-               self.nnearest+=1
-               if self.root is None:
-                       return None, None
-               self.root.count=0
-               node, distsq = self._nearest(self.root, pos, checkempty)
-               self.count+=self.root.count
-               return node,distsq
-               
-       def __str__(self):
-               return str(self.root)
+    """implements a kd-tree"""
+
+    def __init__(self, dim):
+        self.root = None
+        self.nnearest = 0  # number of nearest neighbor queries
+        self.count = 0  # number of nodes visited
+        self.level = 0  # deepest node level
+
+    def resetcounters(self):
+        self.nnearest = 0  # number of nearest neighbor queries
+        self.count = 0  # number of nodes visited
+
+    def _insert(self, node, pos, data):
+        if pos[node.dir] < node.pos[node.dir]:
+            if node.left is None:
+                node.addleft(Node(pos, data))
+                return node.left
+            else:
+                node.rect.extend(pos)
+                return self._insert(node.left, pos, data)
+        else:
+            if node.right is None:
+                node.addright(Node(pos, data))
+                return node.right
+            else:
+                node.rect.extend(pos)
+                return self._insert(node.right, pos, data)
+
+    def insert(self, pos, data):
+        if self.root is None:
+            self.root = Node(pos, data)
+            self.level = self.root.level
+            return self.root
+        else:
+            node = self._insert(self.root, pos, data)
+            if node.level > self.level:
+                self.level = node.level
+            return node
+
+    def _nearest(self, node, pos, checkempty, level=0):
+
+        self.count += 1
+
+        dir = node.dir
+        d = pos[dir] - node.pos[dir]
+
+        result = node
+        distsq = None
+        if checkempty and (node.data is None):
+            result = None
+        else:
+            distsq = node.distance_squared(pos)
+
+        if d <= 0:
+            neartree = node.left
+            fartree = node.right
+        else:
+            neartree = node.right
+            fartree = node.left
+
+        if neartree is not None:
+            nearnode, neardistsq = self._nearest(neartree, pos, checkempty, level + 1)
+            if (result is None) or (neardistsq is not None and neardistsq < distsq):
+                result, distsq = nearnode, neardistsq
+
+        if fartree is not None:
+            if (result is None) or (fartree.rect.distance_squared(pos) < distsq):
+                farnode, fardistsq = self._nearest(fartree, pos, checkempty, level + 1)
+                if (result is None) or (fardistsq is not None and fardistsq < distsq):
+                    result, distsq = farnode, fardistsq
+
+        return result, distsq
+
+    def nearest(self, pos, checkempty=False):
+        self.nnearest += 1
+        if self.root is None:
+            return None, None
+        self.root.count = 0
+        node, distsq = self._nearest(self.root, pos, checkempty)
+        self.count += self.root.count
+        return node, distsq
+
+    def __str__(self):
+        return str(self.root)
 
 if __name__ == "__main__":
 
-       class vector(list):
-
-               def __init__(self, *args):
-                       super().__init__([float(a) for a in args])
-                       
-               def __str__(self):
-                 return "<%.1f %.1f %.1f>"%tuple(self[0:3])
-
-               def __sub__(self,other):
-                 return vector(self[0]-other[0], self[1]-other[1], self[2]-other[2])
-
-               def __add__(self,other):
-                 return vector(self[0]+other[0], self[1]+other[1], self[2]+other[2])
-
-               def __mul__(self,other):
-                 s=sum(self[i]*other[i] for i in (0,1,2))
-                 #print("ds",s,self,other,[self[i]*other[i] for i in (0,1,2)])
-                 return s
-
-               def dot(self,other):
-                 return sum(self[k]*other[k] for k in (0,1,2))
-
-       from random import random,seed, shuffle
-       from time import time
-       import unittest
-
-       class TestVector(unittest.TestCase):
-               def test_ops(self):
-                       v1=vector(1,0,0)
-                       v2=vector(2,1,0)
-                       self.assertAlmostEqual(v1*v2,2.0)
-                       self.assertAlmostEqual(v1.dot(v2),2.0)
-                       v3=vector(-1,-1,0)
-                       self.assertListEqual(v1-v2,v3)
-                       v4=vector(3,1,0)
-                       self.assertListEqual(v1+v2,v4)
-                       
-       class TestHyperRectangle(unittest.TestCase):
-
-               def setUp(self):
-                       self.left=vector(0,0,0)
-                       self.right=vector(1,1,1)
-                       self.left1=vector(-1,0,0)
-                       self.left2=vector(0,-1,0)
-                       self.left3=vector(0,0,-1)
-                       self.right1=vector(2,0,0)
-                       self.right2=vector(0,2,0)
-                       self.right3=vector(0,0,2)
-
-               def test_constructor(self):
-                       hr=Hyperrectangle(3,self.left,self.right)
-                       self.assertListEqual(hr.min,self.left)
-                       self.assertListEqual(hr.max,self.right)
-
-               def test_extend(self):
-                       hr=Hyperrectangle(3,self.left,self.right)
-                       hr.extend(self.left1)
-                       self.assertListEqual(hr.min,[-1,0,0])
-                       self.assertListEqual(hr.max,[1,1,1])
-                       hr.extend(self.left2)
-                       self.assertListEqual(hr.min,[-1,-1,0])
-                       self.assertListEqual(hr.max,[1,1,1])
-                       hr.extend(self.left3)
-                       self.assertListEqual(hr.min,[-1,-1,-1])
-                       self.assertListEqual(hr.max,[1,1,1])
-                       hr.extend(self.right1)
-                       self.assertListEqual(hr.min,[-1,-1,-1])
-                       self.assertListEqual(hr.max,[2,1,1])
-                       hr.extend(self.right2)
-                       self.assertListEqual(hr.min,[-1,-1,-1])
-                       self.assertListEqual(hr.max,[2,2,1])
-                       hr.extend(self.right3)
-                       self.assertListEqual(hr.min,[-1,-1,-1])
-                       self.assertListEqual(hr.max,[2,2,2])
-
-               def test_distance_squared(self):
-                       hr=Hyperrectangle(3,self.left,self.right)
-                       self.assertAlmostEqual(hr.distance_squared(vector(0.5,0.5,0.5)),0.0)
-                       self.assertAlmostEqual(hr.distance_squared(vector(0,0,0)),0.0)
-                       self.assertAlmostEqual(hr.distance_squared(vector(-1,0,0)),1.0)
-                       self.assertAlmostEqual(hr.distance_squared(vector(2,0,0)),1.0)
-                       self.assertAlmostEqual(hr.distance_squared(vector(2,2,2)),3.0)
-                       self.assertAlmostEqual(hr.distance_squared(vector(0.5,2,2)),2.0)
-                       self.assertAlmostEqual(hr.distance_squared(vector(0.5,-1,-1)),2.0)
-                       self.assertAlmostEqual(hr.distance_squared(vector(0.5,0.5,2)),1.0)
-
-       class TestTree(unittest.TestCase):
-
-               def setUp(self):
-                       seed(42)
-                       r=(-1,0,1)
-                       self.points=[vector(k,l,m) for k in r for l in r for m in r]
-                       d=(-0.1,0,0.1)
-                       self.d=[vector(k,l,m) for k in d for l in d for m in d if (k*l*m) != 0]
-                       self.repeats=4
-                       
-               def test_simple(self):
-                       tree=Tree(3)
-                       p1=vector(0,0,0)
-                       p2=vector(-1,0,0)
-                       p3=vector(-1,1,0)
-                       d=vector(0.1,0.1,0.1)
-                       
-                       tree.insert(p1,p1)
-                       node, distsq = tree.nearest(p1)
-                       self.assertListEqual(node.pos,p1)
-                       self.assertAlmostEqual(distsq,0.0)
-                       node, distsq = tree.nearest(p1+d)
-                       self.assertListEqual(node.pos,p1)
-                       self.assertAlmostEqual(distsq,0.03)
-                       
-                       tree.insert(p2,p2)
-                       node, distsq = tree.nearest(p1)
-                       self.assertListEqual(node.pos,p1)
-                       self.assertAlmostEqual(distsq,0.0)
-                       node, distsq = tree.nearest(p1+d)
-                       self.assertListEqual(node.pos,p1)
-                       self.assertAlmostEqual(distsq,0.03)
-                       
-                       node, distsq = tree.nearest(p2)
-                       self.assertListEqual(node.pos,p2)
-                       self.assertAlmostEqual(distsq,0.0)
-                       node, distsq = tree.nearest(p2+d)
-                       self.assertListEqual(node.pos,p2)
-                       self.assertAlmostEqual(distsq,0.03)
-                       
-                       tree.insert(p3,p3)
-                       node, distsq = tree.nearest(p1)
-                       self.assertListEqual(node.pos,p1)
-                       self.assertAlmostEqual(distsq,0.0)
-                       node, distsq = tree.nearest(p1+d)
-                       self.assertListEqual(node.pos,p1)
-                       self.assertAlmostEqual(distsq,0.03)
-                       
-                       node, distsq = tree.nearest(p2)
-                       self.assertListEqual(node.pos,p2)
-                       self.assertAlmostEqual(distsq,0.0)
-                       node, distsq = tree.nearest(p2+d)
-                       self.assertListEqual(node.pos,p2)
-                       self.assertAlmostEqual(distsq,0.03)
-                       
-                       node, distsq = tree.nearest(p3)
-                       self.assertListEqual(node.pos,p3)
-                       self.assertAlmostEqual(distsq,0.0)
-                       node, distsq = tree.nearest(p3+d)
-                       self.assertListEqual(node.pos,p3)
-                       self.assertAlmostEqual(distsq,0.03)
-                       
-               def test_nearest(self):
-                       for n in range(self.repeats):
-                               tree=Tree(3)
-                               shuffle(self.points)
-                               for p in self.points:
-                                       tree.insert(p,p) # data equal to position
-                               
-                               for p in self.points:
-                                       for d in self.d:
-                                               node, distsq = tree.nearest(p+d)
-                                               s="%s %s %s %s\n%s"%(str(p+d), str(p), str(d), str(node.pos),str(tree.root))
-                                               self.assertListEqual(node.pos,p,msg=s)
-                                               self.assertListEqual(node.data,p)
-                                               self.assertAlmostEqual(distsq,d.dot(d),msg=s)
-                               
-                               for p in self.points:
-                                       node, distsq = tree.nearest(p)
-                                       self.assertListEqual(node.pos,p)
-                                       self.assertListEqual(node.data,p)
-                                       self.assertAlmostEqual(distsq,0.0)
-                               
-               def test_nearest_empty(self):
-                       for n in range(self.repeats):
-                               tree=Tree(3)
-                               shuffle(self.points)
-                               for p in self.points:
-                                       tree.insert(p,p) # data equal to position
-                               
-                               for p in self.points:
-                                       for d in self.d:
-                                               node, distsq = tree.nearest(p+d)
-                                               s="%s %s %s %s\n%s"%(str(p+d), str(p), str(d), str(node.pos),str(tree.root))
-                                               self.assertListEqual(node.pos,p,msg=s)
-                                               self.assertListEqual(node.data,p)
-                                               self.assertAlmostEqual(distsq,d.dot(d),msg=s)
-                               
-                               for p in self.points:
-                                       node, distsq = tree.nearest(p)
-                                       self.assertListEqual(node.pos,p)
-                                       self.assertListEqual(node.data,p)
-                                       self.assertAlmostEqual(distsq,0.0)
-
-                               # zeroing out a node should not affect retrieving any other node ...
-                               node , _ = tree.nearest(self.points[-1]) # last point
-                               node.data=None
-                               for p in self.points[:-1]: # all but last
-                                       for d in self.d:
-                                               node, distsq = tree.nearest(p+d)
-                                               s="%s %s %s %s\n%s"%(str(p+d), str(p), str(d), str(node.pos),str(tree.root))
-                                               self.assertListEqual(node.pos,p,msg=s)
-                                               self.assertListEqual(node.data,p)
-                                               self.assertAlmostEqual(distsq,d.dot(d),msg=s)
-                               
-                               for p in self.points[:-1]: # all but last
-                                       node, distsq = tree.nearest(p)
-                                       self.assertListEqual(node.pos,p)
-                                       self.assertListEqual(node.data,p)
-                                       self.assertAlmostEqual(distsq,0.0)
-                               
-                               # ... also, we should be able to retrieve the node anyway ...
-                               node, distsq = tree.nearest(self.points[-1])
-                               self.assertListEqual(node.pos,self.points[-1])
-                               self.assertIsNone(node.data)
-                               self.assertAlmostEqual(distsq,0.0)
-                               
-                               # ... even for points nearby ...
-                               for d in self.d:
-                                       node, distsq = tree.nearest(self.points[-1]+d)
-                                       self.assertEqual(node.pos,self.points[-1])
-                                       self.assertIsNone(node.data)
-                                       self.assertAlmostEqual(distsq,d.dot(d))
-                               
-                               # ... unless we set checkempty
-                               node, distsq = tree.nearest(self.points[-1],checkempty=True)
-                               self.assertNotEqual(node.pos,self.points[-1])
-                               self.assertIsNotNone(node.data)
-                               self.assertAlmostEqual(distsq,1.0) # on a perpendicular grid nearest neighbor is at most 1 away
-                               
-               def test_performance(self):
-                       tree=Tree(3)
-                       tsize=1000
-                       qsize=1000
-                       emptyq=3
-                       
-                       print("<performance test, may take several seconds>")
-                       qpos=[vector(random(),random(),random()) for p in range(qsize)]
-                       for p in range(tsize):
-                               pos=vector(random(),random(),random())
-                               tree.insert(pos,pos)
-                       s=time()
-                       for p in qpos:
-                               node, distsq = tree.nearest(p)
-                       e=time()-s
-                       print("queries|tree size|tree height|empties|query load|query time") 
-                       print("{0:7d}|{2:9d}|{1.level:11d}|      0|{3:10.2f}|{4:10.1f}".format(qsize,tree,tsize,float(tree.count)/qsize,e))
-                       
-                       tree.resetcounters()
-                       empty=[]
-                       for p in range(tsize*9):
-                               pos=vector(random(),random(),random())
-                               tree.insert(pos,pos)
-                               if (p % emptyq ) == 1 :
-                                       empty.append(pos)
-                       s=time()
-                       for p in qpos:
-                               node, distsq = tree.nearest(p)
-                       e2=time()-s
-                       print("{0:7d}|{2:9d}|{1.level:11d}|      0|{3:10.2f}|{4:10.1f}".format(qsize,tree,tsize*10,float(tree.count)/qsize,e2))
-                       
-                       self.assertLess(e2,3*e,msg="a 10x bigger tree shouldn't take more than 3x the time to query")
-
-                       for p in empty:
-                               node, distsq = tree.nearest(p)
-                               node.data = None
-                       tree.resetcounters()
-                       s=time()
-                       for p in qpos:
-                               node, distsq = tree.nearest(p,checkempty=True)
-                       e3=time()-s
-                       print("{0:7d}|{2:9d}|{1.level:11d}|{5:7d}|{3:10.2f}|{4:10.1f}".format(qsize,tree,tsize*10,float(tree.count)/qsize,e3,tsize*10//emptyq))
-                       
-       unittest.main()
\ No newline at end of file
+    class vector(list):
+
+        def __init__(self, *args):
+            super().__init__([float(a) for a in args])
+
+        def __str__(self):
+            return "<%.1f %.1f %.1f>" % tuple(self[0:3])
+
+        def __sub__(self, other):
+            return vector(self[0] - other[0], self[1] - other[1], self[2] - other[2])
+
+        def __add__(self, other):
+            return vector(self[0] + other[0], self[1] + other[1], self[2] + other[2])
+
+        def __mul__(self, other):
+            s = sum(self[i] * other[i] for i in (0, 1, 2))
+            #print("ds",s,self,other,[self[i]*other[i] for i in (0,1,2)])
+            return s
+
+        def dot(self, other):
+            return sum(self[k] * other[k] for k in (0, 1, 2))
+
+    from random import random, seed, shuffle
+    from time import time
+    import unittest
+
+    class TestVector(unittest.TestCase):
+        def test_ops(self):
+            v1 = vector(1, 0, 0)
+            v2 = vector(2, 1, 0)
+            self.assertAlmostEqual(v1 * v2, 2.0)
+            self.assertAlmostEqual(v1.dot(v2), 2.0)
+            v3 = vector(-1, -1, 0)
+            self.assertListEqual(v1 - v2, v3)
+            v4 = vector(3, 1, 0)
+            self.assertListEqual(v1 + v2, v4)
+
+    class TestHyperRectangle(unittest.TestCase):
+
+        def setUp(self):
+            self.left = vector(0, 0, 0)
+            self.right = vector(1, 1, 1)
+            self.left1 = vector(-1, 0, 0)
+            self.left2 = vector(0, -1, 0)
+            self.left3 = vector(0, 0, -1)
+            self.right1 = vector(2, 0, 0)
+            self.right2 = vector(0, 2, 0)
+            self.right3 = vector(0, 0, 2)
+
+        def test_constructor(self):
+            hr = Hyperrectangle(3, self.left, self.right)
+            self.assertListEqual(hr.min, self.left)
+            self.assertListEqual(hr.max, self.right)
+
+        def test_extend(self):
+            hr = Hyperrectangle(3, self.left, self.right)
+            hr.extend(self.left1)
+            self.assertListEqual(hr.min, [-1, 0, 0])
+            self.assertListEqual(hr.max, [1, 1, 1])
+            hr.extend(self.left2)
+            self.assertListEqual(hr.min, [-1, -1, 0])
+            self.assertListEqual(hr.max, [1, 1, 1])
+            hr.extend(self.left3)
+            self.assertListEqual(hr.min, [-1, -1, -1])
+            self.assertListEqual(hr.max, [1, 1, 1])
+            hr.extend(self.right1)
+            self.assertListEqual(hr.min, [-1, -1, -1])
+            self.assertListEqual(hr.max, [2, 1, 1])
+            hr.extend(self.right2)
+            self.assertListEqual(hr.min, [-1, -1, -1])
+            self.assertListEqual(hr.max, [2, 2, 1])
+            hr.extend(self.right3)
+            self.assertListEqual(hr.min, [-1, -1, -1])
+            self.assertListEqual(hr.max, [2, 2, 2])
+
+        def test_distance_squared(self):
+            hr = Hyperrectangle(3, self.left, self.right)
+            self.assertAlmostEqual(hr.distance_squared(vector(0.5, 0.5, 0.5)), 0.0)
+            self.assertAlmostEqual(hr.distance_squared(vector(0, 0, 0)), 0.0)
+            self.assertAlmostEqual(hr.distance_squared(vector(-1, 0, 0)), 1.0)
+            self.assertAlmostEqual(hr.distance_squared(vector(2, 0, 0)), 1.0)
+            self.assertAlmostEqual(hr.distance_squared(vector(2, 2, 2)), 3.0)
+            self.assertAlmostEqual(hr.distance_squared(vector(0.5, 2, 2)), 2.0)
+            self.assertAlmostEqual(hr.distance_squared(vector(0.5, -1, -1)), 2.0)
+            self.assertAlmostEqual(hr.distance_squared(vector(0.5, 0.5, 2)), 1.0)
+
+    class TestTree(unittest.TestCase):
+
+        def setUp(self):
+            seed(42)
+            r = (-1, 0, 1)
+            self.points = [vector(k, l, m) for k in r for l in r for m in r]
+            d = (-0.1, 0, 0.1)
+            self.d = [vector(k, l, m) for k in d for l in d for m in d if (k * l * m) != 0]
+            self.repeats = 4
+
+        def test_simple(self):
+            tree = Tree(3)
+            p1 = vector(0, 0, 0)
+            p2 = vector(-1, 0, 0)
+            p3 = vector(-1, 1, 0)
+            d = vector(0.1, 0.1, 0.1)
+
+            tree.insert(p1, p1)
+            node, distsq = tree.nearest(p1)
+            self.assertListEqual(node.pos, p1)
+            self.assertAlmostEqual(distsq, 0.0)
+            node, distsq = tree.nearest(p1 + d)
+            self.assertListEqual(node.pos, p1)
+            self.assertAlmostEqual(distsq, 0.03)
+
+            tree.insert(p2, p2)
+            node, distsq = tree.nearest(p1)
+            self.assertListEqual(node.pos, p1)
+            self.assertAlmostEqual(distsq, 0.0)
+            node, distsq = tree.nearest(p1 + d)
+            self.assertListEqual(node.pos, p1)
+            self.assertAlmostEqual(distsq, 0.03)
+
+            node, distsq = tree.nearest(p2)
+            self.assertListEqual(node.pos, p2)
+            self.assertAlmostEqual(distsq, 0.0)
+            node, distsq = tree.nearest(p2 + d)
+            self.assertListEqual(node.pos, p2)
+            self.assertAlmostEqual(distsq, 0.03)
+
+            tree.insert(p3, p3)
+            node, distsq = tree.nearest(p1)
+            self.assertListEqual(node.pos, p1)
+            self.assertAlmostEqual(distsq, 0.0)
+            node, distsq = tree.nearest(p1 + d)
+            self.assertListEqual(node.pos, p1)
+            self.assertAlmostEqual(distsq, 0.03)
+
+            node, distsq = tree.nearest(p2)
+            self.assertListEqual(node.pos, p2)
+            self.assertAlmostEqual(distsq, 0.0)
+            node, distsq = tree.nearest(p2 + d)
+            self.assertListEqual(node.pos, p2)
+            self.assertAlmostEqual(distsq, 0.03)
+
+            node, distsq = tree.nearest(p3)
+            self.assertListEqual(node.pos, p3)
+            self.assertAlmostEqual(distsq, 0.0)
+            node, distsq = tree.nearest(p3 + d)
+            self.assertListEqual(node.pos, p3)
+            self.assertAlmostEqual(distsq, 0.03)
+
+        def test_nearest(self):
+            for n in range(self.repeats):
+                tree = Tree(3)
+                shuffle(self.points)
+                for p in self.points:
+                    tree.insert(p, p)  # data equal to position
+
+                for p in self.points:
+                    for d in self.d:
+                        node, distsq = tree.nearest(p + d)
+                        s = "%s %s %s %s\n%s" % (str(p + d), str(p), str(d), str(node.pos), str(tree.root))
+                        self.assertListEqual(node.pos, p, msg=s)
+                        self.assertListEqual(node.data, p)
+                        self.assertAlmostEqual(distsq, d.dot(d), msg=s)
+
+                for p in self.points:
+                    node, distsq = tree.nearest(p)
+                    self.assertListEqual(node.pos, p)
+                    self.assertListEqual(node.data, p)
+                    self.assertAlmostEqual(distsq, 0.0)
+
+        def test_nearest_empty(self):
+            for n in range(self.repeats):
+                tree = Tree(3)
+                shuffle(self.points)
+                for p in self.points:
+                    tree.insert(p, p)  # data equal to position
+
+                for p in self.points:
+                    for d in self.d:
+                        node, distsq = tree.nearest(p + d)
+                        s = "%s %s %s %s\n%s" % (str(p + d), str(p), str(d), str(node.pos), str(tree.root))
+                        self.assertListEqual(node.pos, p, msg=s)
+                        self.assertListEqual(node.data, p)
+                        self.assertAlmostEqual(distsq, d.dot(d), msg=s)
+
+                for p in self.points:
+                    node, distsq = tree.nearest(p)
+                    self.assertListEqual(node.pos, p)
+                    self.assertListEqual(node.data, p)
+                    self.assertAlmostEqual(distsq, 0.0)
+
+                # zeroing out a node should not affect retrieving any other node ...
+                node, _ = tree.nearest(self.points[-1])  # last point
+                node.data = None
+                for p in self.points[:-1]:  # all but last
+                    for d in self.d:
+                        node, distsq = tree.nearest(p + d)
+                        s = "%s %s %s %s\n%s" % (str(p + d), str(p), str(d), str(node.pos), str(tree.root))
+                        self.assertListEqual(node.pos, p, msg=s)
+                        self.assertListEqual(node.data, p)
+                        self.assertAlmostEqual(distsq, d.dot(d), msg=s)
+
+                for p in self.points[:-1]:  # all but last
+                    node, distsq = tree.nearest(p)
+                    self.assertListEqual(node.pos, p)
+                    self.assertListEqual(node.data, p)
+                    self.assertAlmostEqual(distsq, 0.0)
+
+                # ... also, we should be able to retrieve the node anyway ...
+                node, distsq = tree.nearest(self.points[-1])
+                self.assertListEqual(node.pos, self.points[-1])
+                self.assertIsNone(node.data)
+                self.assertAlmostEqual(distsq, 0.0)
+
+                # ... even for points nearby ...
+                for d in self.d:
+                    node, distsq = tree.nearest(self.points[-1] + d)
+                    self.assertEqual(node.pos, self.points[-1])
+                    self.assertIsNone(node.data)
+                    self.assertAlmostEqual(distsq, d.dot(d))
+
+                # ... unless we set checkempty
+                node, distsq = tree.nearest(self.points[-1], checkempty=True)
+                self.assertNotEqual(node.pos, self.points[-1])
+                self.assertIsNotNone(node.data)
+                self.assertAlmostEqual(distsq, 1.0)  # on a perpendicular grid nearest neighbor is at most 1 away
+
+        def test_performance(self):
+            tree = Tree(3)
+            tsize = 1000
+            qsize = 1000
+            emptyq = 3
+
+            print("<performance test, may take several seconds>")
+            qpos = [vector(random(), random(), random()) for p in range(qsize)]
+            for p in range(tsize):
+                pos = vector(random(), random(), random())
+                tree.insert(pos, pos)
+            s = time()
+            for p in qpos:
+                node, distsq = tree.nearest(p)
+            e = time() - s
+            print("queries|tree size|tree height|empties|query load|query time")
+            print("{0:7d}|{2:9d}|{1.level:11d}|      0|{3:10.2f}|{4:10.1f}".format(qsize, tree, tsize, float(tree.count) / qsize, e))
+
+            tree.resetcounters()
+            empty = []
+            for p in range(tsize * 9):
+                pos = vector(random(), random(), random())
+                tree.insert(pos, pos)
+                if (p % emptyq) == 1:
+                    empty.append(pos)
+            s = time()
+            for p in qpos:
+                node, distsq = tree.nearest(p)
+            e2 = time() - s
+            print("{0:7d}|{2:9d}|{1.level:11d}|      0|{3:10.2f}|{4:10.1f}".format(qsize, tree, tsize * 10, float(tree.count) / qsize, e2))
+
+            self.assertLess(e2, 3 * e, msg="a 10x bigger tree shouldn't take more than 3x the time to query")
+
+            for p in empty:
+                node, distsq = tree.nearest(p)
+                node.data = None
+            tree.resetcounters()
+            s = time()
+            for p in qpos:
+                node, distsq = tree.nearest(p, checkempty=True)
+            e3 = time() - s
+            print("{0:7d}|{2:9d}|{1.level:11d}|{5:7d}|{3:10.2f}|{4:10.1f}".format(qsize, tree, tsize * 10, float(tree.count) / qsize, e3, tsize * 10 // emptyq))
+
+    unittest.main()
index 55cac59..f9d2c6b 100644 (file)
 #
 # ##### END GPL LICENSE BLOCK #####
 
-
+# <pep8 compliant>
 
 from collections import defaultdict as dd
-from random import random,seed,expovariate
-from math import sqrt,pow,sin,cos
+from random import random, seed, expovariate
+from math import sqrt, pow, sin, cos
 from functools import partial
 
 from mathutils import Vector
 
 from .kdtree import Tree
 
+
 class Branchpoint:
-       def __init__(self, p, parent):
-               self.v=Vector(p)
-               self.parent = parent
-               self.connections = 0
-               self.apex = None
-               self.shoot = None
-               self.index = None
-
-def sphere(r,p):
-       r2 = r*r
-       while True:
-               x = (random()*2-1)*r
-               y = (random()*2-1)*r
-               z = (random()*2-1)*r
-               if x*x+y*y+z*z <= r2:
-                       yield p+Vector((x,y,z))
+    def __init__(self, p, parent):
+        self.v = Vector(p)
+        self.parent = parent
+        self.connections = 0
+        self.apex = None
+        self.shoot = None
+        self.index = None
+
+
+def sphere(r, p):
+    r2 = r * r
+    while True:
+        x = (random() * 2 - 1) * r
+        y = (random() * 2 - 1) * r
+        z = (random() * 2 - 1) * r
+        if x * x + y * y + z * z <= r2:
+            yield p + Vector((x, y, z))
+
 
 class SCA:
-       def __init__(self,NENDPOINTS = 100,d = 0.3,NBP = 2000, KILLDIST = 5, INFLUENCE = 15, SEED=42, volume=partial(sphere,5,Vector((0,0,8))), TROPISM=0.0, exclude=lambda p: False,
-               startingpoints=[]):
-               seed(SEED)
-               self.d = d
-               self.NBP = NBP
-               self.KILLDIST = KILLDIST*KILLDIST*d*d
-               self.INFLUENCE = INFLUENCE*INFLUENCE*d*d
-               self.TROPISM = TROPISM
-               self.exclude = exclude
-               self.endpoints = []
-
-               self.volumepoint=volume()
-               for i in range(NENDPOINTS):
-                       self.endpoints.append(next(self.volumepoint))
-       
-               self.branchpoints = [ Branchpoint((0,0,0),None) ] if len(startingpoints)==0 else startingpoints
-               
-       def iterate(self, newendpointsper1000=0, maxtime=0.0): # maxtime still ignored for now
-               
-               endpointsadded=0.0
-               niterations=0.0
-               newendpointsper1000 /= 1000.0
-               t=expovariate(newendpointsper1000) if newendpointsper1000 > 0.0 else 1 # time to the first new 'endpoint add event'
-               
-               while self.NBP>0 and (len(self.endpoints)>0):
-                       self.NBP -= 1
-                       closestsendpoints=dd(list)
-
-                       kill = set()
-                       
-                       for ei,e in enumerate(self.endpoints):
-                               distance = None
-                               closestbp = None
-                               for bi,b in enumerate(self.branchpoints):
-                                       ddd = b.v-e
-                                       ddd = ddd.dot(ddd)
-                                       if ddd < self.KILLDIST:
-                                               kill.add(ei)
-                                       elif (ddd<self.INFLUENCE and b.shoot is None) and ((distance is None) or (ddd < distance)):
-                                               closestbp = bi
-                                               distance = ddd
-                               if not (closestbp is None):
-                                       closestsendpoints[closestbp].append(ei)
-                       
-                       if len(closestsendpoints)<1:
-                               break
-                       
-                       for bi in closestsendpoints:
-                               sd=Vector((0,0,0))
-                               n=0
-                               for ei in closestsendpoints[bi]:
-                                       dv=self.branchpoints[bi].v-self.endpoints[ei]
-                                       ll=sqrt(dv.dot(dv))
-                                       sd-=dv/ll
-                                       n+=1
-                               sd/=n
-                               ll=sqrt(sd.dot(sd))
-                               # don't know if this assumption is true:
-                               # if the unnormalised direction is very small, the endpoints are nearly coplanar/colinear and at roughly the same distance
-                               # so no endpoints will be killed and we might end up adding the same branch again and again
-                               if ll < 1e-3 :
-                                       #print('SD very small')
-                                       continue
-                                       
-                               sd/=ll
-                               sd[2]+=self.TROPISM
-                               ll=sqrt(sd.dot(sd))
-                               sd/=ll
-                       
-                               newp = self.branchpoints[bi].v+sd*self.d
-                               # the assumption we made earlier is not suffucient to prevent adding the same branch so we need an expensive check:
-                               tooclose = False
-                               for dbi in self.branchpoints:
-                                       dddd = newp - dbi.v
-                                       if dddd.dot(dddd) < 1e-3 :
-                                               #print('BP to close to another')
-                                               tooclose = True
-                               if tooclose : continue
-                               
-                               if not self.exclude(newp):
-                                       bp = Branchpoint(newp,bi)
-                                       self.branchpoints.append(bp)
-                                       nbpi = len(self.branchpoints)-1
-                                       bp = self.branchpoints[bi]
-                                       bp.connections+=1
-                                       if bp.apex is None:
-                                               bp.apex = nbpi
-                                       else:
-                                               bp.shoot = nbpi
-                                       while not (bp.parent is None):
-                                               bp = self.branchpoints[bp.parent]
-                                               bp.connections+=1
-                                       
-                       self.endpoints = [ep for ei,ep in enumerate(self.endpoints) if not(ei in kill)]
-               
-                       if newendpointsper1000 > 0.0:
-                               # generate new endpoints with a poisson process
-                               # when we first arrive here, t already hold the time to the first event
-                               niterations+=1
-                               while t < niterations: # we keep on adding endpoints as long as the next event still happens within this iteration
-                                       self.endpoints.append(next(self.volumepoint))
-                                       endpointsadded+=1
-                                       t+=expovariate(newendpointsper1000) # time to new 'endpoint add event'
-                               
-               #if newendpointsper1000 > 0.0:
-                       #print("newendpoints/iteration %.3f, actual %.3f in %5.f iterations"%(newendpointsper1000,endpointsadded/niterations,niterations))
-
-       def iterate2(self, newendpointsper1000=0, maxtime=0.0): # maxtime still ignored for now
-               """iterate using a kdtree fr the branchpoints"""
-               endpointsadded=0.0
-               niterations=0.0
-               newendpointsper1000 /= 1000.0
-               t=expovariate(newendpointsper1000) if newendpointsper1000 > 0.0 else 1 # time to the first new 'endpoint add event'
-               
-               tree=Tree(3)
-               for bpi,bp in enumerate(self.branchpoints):
-                       bp.index=bpi
-                       tree.insert(bp.v,bp)
-
-               while self.NBP>0 and (len(self.endpoints)>0):
-                       self.NBP -= 1
-                       closestsendpoints=dd(list)
-
-                       kill = set()
-                       
-                       for ei,e in enumerate(self.endpoints):
-                               distance = None
-                               closestbp, distance = tree.nearest(e, checkempty=True) # ignore forks, they have .data set to None
-                               if closestbp is not None:
-                                       if distance < self.KILLDIST:
-                                               kill.add(ei)
-                                       elif distance<self.INFLUENCE : 
-                                               closestsendpoints[closestbp.data.index].append(ei)
-                       
-                       if len(closestsendpoints)<1:
-                               break
-                       
-                       for bi in closestsendpoints:
-                               sd=Vector((0,0,0))
-                               n=0
-                               for ei in closestsendpoints[bi]:
-                                       dv=self.branchpoints[bi].v-self.endpoints[ei]
-                                       dv.normalize()
-                                       sd-=dv
-                                       n+=1
-                               sd/=n
-                               ll=sd.length_squared
-                               # don't know if this assumption is true:
-                               # if the unnormalised direction is very small, the endpoints are nearly coplanar/colinear and at roughly the same distance
-                               # so no endpoints will be killed and we might end up adding the same branch again and again
-                               if ll < 1e-3 :
-                                       #print('SD very small')
-                                       continue
-                                       
-                               sd/=ll
-                               sd[2]+=self.TROPISM
-                               sd.normalize()
-                       
-                               newp = self.branchpoints[bi].v+sd*self.d
-                               # the assumption we made earlier is not suffucient to prevent adding the same branch so we need an expensive check:
-                               _, dddd = tree.nearest(newp) # no checkempty here, we want to check for forks as well
-                               if dddd < 1e-3 :
-                                               #print('BP to close to another')
-                                               continue
-                               
-                               if not self.exclude(newp):
-                                       bp = Branchpoint(newp,bi)
-                                       self.branchpoints.append(bp)
-                                       nbpi = len(self.branchpoints)-1
-                                       bp.index=nbpi
-                                       tree.insert(bp.v,bp)
-                                       bp = self.branchpoints[bi]
-                                       bp.connections+=1
-                                       if bp.apex is None:
-                                               bp.apex = nbpi
-                                       else:
-                                               bp.shoot = nbpi
-                                               node,_ = tree.nearest(bp.v)
-                                               node.data = None # signal that this is a fork so that we might ignore it when searching for nearest neighbors
-                                       while not (bp.parent is None):
-                                               bp = self.branchpoints[bp.parent]
-                                               bp.connections+=1
-                                       
-                       self.endpoints = [ep for ei,ep in enumerate(self.endpoints) if not(ei in kill)]
-               
-                       if newendpointsper1000 > 0.0:
-                               # generate new endpoints with a poisson process
-                               # when we first arrive here, t already hold the time to the first event
-                               niterations+=1
-                               while t < niterations: # we keep on adding endpoints as long as the next event still happens within this iteration
-                                       self.endpoints.append(next(self.volumepoint))
-                                       endpointsadded+=1
-                                       t+=expovariate(newendpointsper1000) # time to new 'endpoint add event'
-                               
\ No newline at end of file
+    def __init__(self, NENDPOINTS=100, d=0.3, NBP=2000, KILLDIST=5, INFLUENCE=15, SEED=42, volume=partial(sphere, 5, Vector((0, 0, 8))), TROPISM=0.0, exclude=lambda p: False,
+        startingpoints=[]):
+        seed(SEED)
+        self.d = d
+        self.NBP = NBP
+        self.KILLDIST = KILLDIST * KILLDIST * d * d
+        self.INFLUENCE = INFLUENCE * INFLUENCE * d * d
+        self.TROPISM = TROPISM
+        self.exclude = exclude
+        self.endpoints = []
+
+        self.volumepoint = volume()
+        for i in range(NENDPOINTS):
+            self.endpoints.append(next(self.volumepoint))
+
+        self.branchpoints = [Branchpoint((0, 0, 0), None)] if len(startingpoints) == 0 else startingpoints
+
+    def iterate(self, newendpointsper1000=0, maxtime=0.0):  # maxtime still ignored for now
+
+        endpointsadded = 0.0
+        niterations = 0.0
+        newendpointsper1000 /= 1000.0
+        t = expovariate(newendpointsper1000) if newendpointsper1000 > 0.0 else 1  # time to the first new 'endpoint add event'
+
+        while self.NBP > 0 and (len(self.endpoints) > 0):
+            self.NBP -= 1
+            closestsendpoints = dd(list)
+
+            kill = set()
+
+            for ei, e in enumerate(self.endpoints):
+                distance = None
+                closestbp = None
+                for bi, b in enumerate(self.branchpoints):
+                    ddd = b.v - e
+                    ddd = ddd.dot(ddd)
+                    if ddd < self.KILLDIST:
+                        kill.add(ei)
+                    elif (ddd < self.INFLUENCE and b.shoot is None) and ((distance is None) or (ddd < distance)):
+                        closestbp = bi
+                        distance = ddd
+                if not (closestbp is None):
+                    closestsendpoints[closestbp].append(ei)
+
+            if len(closestsendpoints) < 1:
+                break
+
+            for bi in closestsendpoints:
+                sd = Vector((0, 0, 0))
+                n = 0
+                for ei in closestsendpoints[bi]:
+                    dv = self.branchpoints[bi].v - self.endpoints[ei]
+                    ll = sqrt(dv.dot(dv))
+                    sd -= dv / ll
+                    n += 1
+                sd /= n
+                ll = sqrt(sd.dot(sd))
+                # don't know if this assumption is true:
+                # if the unnormalised direction is very small, the endpoints are nearly coplanar/colinear and at roughly the same distance
+                # so no endpoints will be killed and we might end up adding the same branch again and again
+                if ll < 1e-3:
+                    #print('SD very small')
+                    continue
+
+                sd /= ll
+                sd[2] += self.TROPISM
+                ll = sqrt(sd.dot(sd))
+                sd /= ll
+
+                newp = self.branchpoints[bi].v + sd * self.d
+                # the assumption we made earlier is not suffucient to prevent adding the same branch so we need an expensive check:
+                tooclose = False
+                for dbi in self.branchpoints:
+                    dddd = newp - dbi.v
+                    if dddd.dot(dddd) < 1e-3:
+                        #print('BP to close to another')
+                        tooclose = True
+                if tooclose:
+                    continue
+
+                if not self.exclude(newp):
+                    bp = Branchpoint(newp, bi)
+                    self.branchpoints.append(bp)
+                    nbpi = len(self.branchpoints) - 1
+                    bp = self.branchpoints[bi]
+                    bp.connections += 1
+                    if bp.apex is None:
+                        bp.apex = nbpi
+                    else:
+                        bp.shoot = nbpi
+                    while not (bp.parent is None):
+                        bp = self.branchpoints[bp.parent]
+                        bp.connections += 1
+
+            self.endpoints = [ep for ei, ep in enumerate(self.endpoints) if not(ei in kill)]
+
+            if newendpointsper1000 > 0.0:
+                # generate new endpoints with a poisson process
+                # when we first arrive here, t already hold the time to the first event
+                niterations += 1
+                while t < niterations:  # we keep on adding endpoints as long as the next event still happens within this iteration
+                    self.endpoints.append(next(self.volumepoint))
+                    endpointsadded += 1
+                    t += expovariate(newendpointsper1000)  # time to new 'endpoint add event'
+
+        #if newendpointsper1000 > 0.0:
+            #print("newendpoints/iteration %.3f, actual %.3f in %5.f iterations"%(newendpointsper1000,endpointsadded/niterations,niterations))
+
+    def iterate2(self, newendpointsper1000=0, maxtime=0.0):  # maxtime still ignored for now
+        """iterate using a kdtree fr the branchpoints"""
+        endpointsadded = 0.0
+        niterations = 0.0
+        newendpointsper1000 /= 1000.0
+        t = expovariate(newendpointsper1000) if newendpointsper1000 > 0.0 else 1  # time to the first new 'endpoint add event'
+
+        tree = Tree(3)
+        for bpi, bp in enumerate(self.branchpoints):
+            bp.index = bpi
+            tree.insert(bp.v, bp)
+
+        while self.NBP > 0 and (len(self.endpoints) > 0):
+            self.NBP -= 1
+            closestsendpoints = dd(list)
+
+            kill = set()
+
+            for ei, e in enumerate(self.endpoints):
+                distance = None
+                closestbp, distance = tree.nearest(e, checkempty=True)  # ignore forks, they have .data set to None
+                if closestbp is not None:
+                    if distance < self.KILLDIST:
+                        kill.add(ei)
+                    elif distance < self.INFLUENCE:
+                        closestsendpoints[closestbp.data.index].append(ei)
+
+            if len(closestsendpoints) < 1:
+                break
+
+            for bi in closestsendpoints:
+                sd = Vector((0, 0, 0))
+                n = 0
+                for ei in closestsendpoints[bi]:
+                    dv = self.branchpoints[bi].v - self.endpoints[ei]
+                    dv.normalize()
+                    sd -= dv
+                    n += 1
+                sd /= n
+                ll = sd.length_squared
+                # don't know if this assumption is true:
+                # if the unnormalised direction is very small, the endpoints are nearly coplanar/colinear and at roughly the same distance
+                # so no endpoints will be killed and we might end up adding the same branch again and again
+                if ll < 1e-3:
+                    #print('SD very small')
+                    continue
+
+                sd /= ll
+                sd[2] += self.TROPISM
+                sd.normalize()
+
+                newp = self.branchpoints[bi].v + sd * self.d
+                # the assumption we made earlier is not suffucient to prevent adding the same branch so we need an expensive check:
+                _, dddd = tree.nearest(newp)  # no checkempty here, we want to check for forks as well
+                if dddd < 1e-3:
+                        #print('BP to close to another')
+                        continue
+
+                if not self.exclude(newp):
+                    bp = Branchpoint(newp, bi)
+                    self.branchpoints.append(bp)
+                    nbpi = len(self.branchpoints) - 1
+                    bp.index = nbpi
+                    tree.insert(bp.v, bp)
+                    bp = self.branchpoints[bi]
+                    bp.connections += 1
+                    if bp.apex is None:
+                        bp.apex = nbpi
+                    else:
+                        bp.shoot = nbpi
+                        node, _ = tree.nearest(bp.v)
+                        node.data = None  # signal that this is a fork so that we might ignore it when searching for nearest neighbors
+                    while not (bp.parent is None):
+                        bp = self.branchpoints[bp.parent]
+                        bp.connections += 1
+
+            self.endpoints = [ep for ei, ep in enumerate(self.endpoints) if not(ei in kill)]
+
+            if newendpointsper1000 > 0.0:
+                # generate new endpoints with a poisson process
+                # when we first arrive here, t already hold the time to the first event
+                niterations += 1
+                while t < niterations:  # we keep on adding endpoints as long as the next event still happens within this iteration
+                    self.endpoints.append(next(self.volumepoint))
+                    endpointsadded += 1
+                    t += expovariate(newendpointsper1000)  # time to new 'endpoint add event'
index d704d4a..98bea95 100644 (file)
@@ -1,12 +1,12 @@
 # ##### BEGIN GPL LICENSE BLOCK #####
 #
-#  SCA Tree Generator, a Blender addon
+#  SCA Tree Generator,  a Blender addon
 #  (c) 2013 Michel J. Anders (varkenvarken)
 #
 #  This program is free software; you can redistribute it and/or
 #  modify it under the terms of the GNU General Public License
 #  as published by the Free Software Foundation; either version 2
-#  of the License, or (at your option) any later version.
+#  of the License,  or (at your option) any later version.
 #
 #  This program is distributed in the hope that it will be useful,
 #  but WITHOUT ANY WARRANTY; without even the implied warranty of
 #  GNU General Public License for more details.
 #
 #  You should have received a copy of the GNU General Public License
-#  along with this program; if not, write to the Free Software Foundation,
-#  Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+#  along with this program; if not,  write to the Free Software Foundation,
+#  Inc.,  51 Franklin Street,  Fifth Floor,  Boston,  MA 02110-1301,  USA.
 #
 # ##### END GPL LICENSE BLOCK #####
 
+# <pep8 compliant>
+
 from math import pi
 from mathutils import Quaternion
 
-rot120 = 2*pi/3
+rot120 = 2 * pi / 3
+
 
 def rot(point, axis, angle):
-       q=Quaternion(axis,angle)
-       P=point.copy()
-       P.rotate(q)
-       #print(point,P)
-       return P
-
-def vertexnormal(d1,d2,d3):
-  n1=d1.cross(d2).normalized()
-  n2=d2.cross(d3).normalized()
-  n3=d3.cross(d1).normalized()
-  n=(n1+n2+n3).normalized()
-  if (d1+d2+d3).dot(n)>0:
-    return -n
-  return n
-  
-def simplefork2(p0,p1,p2,p3,r0,r1,r2,r3):
-       d1=p1-p0
-       d2=p2-p0
-       d3=p3-p0
-       #print(d1,d2,d3)
-       n=vertexnormal(d1,d2,d3)
-       #print(n)
-
-       pp1=p0+d1/3
-       n1a=r1*n
-       n1b=rot(n1a,d1, rot120)
-       n1c=rot(n1a,d1,-rot120)
-       v1a=pp1+n1a
-       v1b=pp1+n1b
-       v1c=pp1+n1c
-
-       pp2=p0+d2/3
-       n2a=r2*n
-       n2b=rot(n2a,d2, rot120)
-       n2c=rot(n2a,d2,-rot120)
-       v2a=pp2+n2a
-       v2b=pp2+n2b
-       v2c=pp2+n2c
-
-       pp3=p0+d3/3
-       n3a=r3*n
-       n3b=rot(n3a,d3, rot120)
-       n3c=rot(n3a,d3,-rot120)
-       v3a=pp3+n3a
-       v3b=pp3+n3b
-       v3c=pp3+n3c
-
-       n0a=n*r0
-       v0a=p0+n0a
-       v0c=p0-d3.normalized()*r0-n0a/3
-       v0d=p0-d1.normalized()*r0-n0a/3
-       v0b=p0-d2.normalized()*r0-n0a/3
-
-       #v0b=p0+(n1b+n2c)/2
-       #v0d=p0+(n2b+n3c)/2
-       #v0c=p0+(n3b+n1c)/2
-
-       verts=(v1a,v1b,v1c,v2a,v2b,v2c,v3a,v3b,v3c,v0a,v0b,v0c,v0d)
-       faces=( (0,1,10,9),(1,2,11,10),(2,0,9,11),#chck
-               (3,4,11,9),(4,5,12,11),(5,3,9,12),#chck
-               
-               (6,7,12,9),
-               (7,8,10,12),
-               (8,6,9,10),
-               
-               (10,11,12))
-
-       return verts,faces
-
-def simplefork(p0,p1,p2,p3,r0,r1,r2,r3):
-       d1=p1-p0
-       d2=p2-p0
-       d3=p3-p0
-       #print(d1,d2,d3)
-       n=-vertexnormal(d1,d2,d3)
-       #print(n)
-
-       # the central tetrahedron
-       n0a=n*r0*0.3
-       v0a=n0a
-       v0b=-d1/6-n0a/2
-       v0c=-d2/6-n0a/2
-       v0d=-d3/6-n0a/2
-
-       n1=v0a+v0c+v0d
-       n2=v0a+v0b+v0d
-       n3=v0a+v0b+v0c
-
-       q1=n1.rotation_difference(d1)  
-       q2=n2.rotation_difference(d2)  
-       q3=n3.rotation_difference(d3)  
-
-       pp1=p0+d1/3
-       v1a=v0a.copy()
-       v1b=v0c.copy()
-       v1c=v0d.copy()
-       v1a.rotate(q1)
-       v1b.rotate(q1)
-       v1c.rotate(q1)
-       v1a+=pp1
-       v1b+=pp1
-       v1c+=pp1
-
-       pp2=p0+d2/3
-       v2a=v0a.copy()
-       v2b=v0b.copy()
-       v2c=v0d.copy()
-       v2a.rotate(q2)
-       v2b.rotate(q2)
-       v2c.rotate(q2)
-       v2a+=pp2
-       v2b+=pp2
-       v2c+=pp2
-
-       pp3=p0+d3/3
-       v3a=v0a.copy()
-       v3b=v0b.copy()
-       v3c=v0c.copy()
-       v3a.rotate(q3)
-       v3b.rotate(q3)
-       v3c.rotate(q3)
-       v3a+=pp3
-       v3b+=pp3
-       v3c+=pp3
-
-       v0a += p0
-       v0b += p0
-       v0c += p0
-       v0d += p0
-       
-       verts=(v1a,v1b,v1c,v2a,v2b,v2c,v3a,v3b,v3c,v0a,v0b,v0c,v0d)
-       faces=( 
-               #(1,2,12,11),
-               #(9,12,2,0),
-               #(11,9,0,1),
-               
-               #(5,4,10,12),
-               #(4,3,9,10),
-               #(3,5,12,9),
-               
-               (8,7,11,10),
-               (7,5,9,11),
-               (6,8,10,9),
-               
-               (10,11,12))
-
-       return verts,faces
+    q = Quaternion(axis, angle)
+    P = point.copy()
+    P.rotate(q)
+    #print(point, P)
+    return P
+
+
+def vertexnormal(d1, d2, d3):
+    n1 = d1.cross(d2).normalized()
+    n2 = d2.cross(d3).normalized()
+    n3 = d3.cross(d1).normalized()
+    n = (n1 + n2 + n3).normalized()
+    if (d1 + d2 + d3).dot(n) > 0:
+        return -n
+    return n
+
+
+def simplefork2(p0, p1, p2, p3, r0, r1, r2, r3):
+    d1 = p1 - p0
+    d2 = p2 - p0
+    d3 = p3 - p0
+    #print(d1, d2, d3)
+    n = vertexnormal(d1, d2, d3)
+    #print(n)
+
+    pp1 = p0 + d1 / 3
+    n1a = r1 * n
+    n1b = rot(n1a, d1, rot120)
+    n1c = rot(n1a, d1, -rot120)
+    v1a = pp1 + n1a
+    v1b = pp1 + n1b
+    v1c = pp1 + n1c
+
+    pp2 = p0 + d2 / 3
+    n2a = r2 * n
+    n2b = rot(n2a, d2, rot120)
+    n2c = rot(n2a, d2, -rot120)
+    v2a = pp2 + n2a
+    v2b = pp2 + n2b
+    v2c = pp2 + n2c
+
+    pp3 = p0 + d3 / 3
+    n3a = r3 * n
+    n3b = rot(n3a, d3, rot120)
+    n3c = rot(n3a, d3, -rot120)
+    v3a = pp3 + n3a
+    v3b = pp3 + n3b
+    v3c = pp3 + n3c
+
+    n0a = n * r0
+    v0a = p0 + n0a
+    v0c = p0 - d3.normalized() * r0 - n0a / 3
+    v0d = p0 - d1.normalized() * r0 - n0a / 3
+    v0b = p0 - d2.normalized() * r0 - n0a / 3
+
+    #v0b=p0+(n1b+n2c)/2
+    #v0d=p0+(n2b+n3c)/2
+    #v0c=p0+(n3b+n1c)/2
+
+    verts = (v1a, v1b, v1c, v2a, v2b, v2c, v3a, v3b, v3c, v0a, v0b, v0c, v0d)
+    faces = ((0, 1, 10, 9), (1, 2, 11, 10), (2, 0, 9, 11),  # chck
+        (3, 4, 11, 9), (4, 5, 12, 11), (5, 3, 9, 12),  # chck
+
+        (6, 7, 12, 9),
+        (7, 8, 10, 12),
+        (8, 6, 9, 10),
+
+        (10, 11, 12))
+
+    return verts, faces
+
+
+def simplefork(p0, p1, p2, p3, r0, r1, r2, r3):
+    d1 = p1 - p0
+    d2 = p2 - p0
+    d3 = p3 - p0
+    #print(d1, d2, d3)
+    n = -vertexnormal(d1, d2, d3)
+    #print(n)
+
+    # the central tetrahedron
+    n0a = n * r0 * 0.3
+    v0a = n0a
+    v0b = -d1 / 6 - n0a / 2
+    v0c = -d2 / 6 - n0a / 2
+    v0d = -d3 / 6 - n0a / 2
+
+    n1 = v0a + v0c + v0d
+    n2 = v0a + v0b + v0d
+    n3 = v0a + v0b + v0c
+
+    q1 = n1.rotation_difference(d1)
+    q2 = n2.rotation_difference(d2)
+    q3 = n3.rotation_difference(d3)
+
+    pp1 = p0 + d1 / 3
+    v1a = v0a.copy()
+    v1b = v0c.copy()
+    v1c = v0d.copy()
+    v1a.rotate(q1)
+    v1b.rotate(q1)
+    v1c.rotate(q1)
+    v1a += pp1
+    v1b += pp1
+    v1c += pp1
+
+    pp2 = p0 + d2 / 3
+    v2a = v0a.copy()
+    v2b = v0b.copy()
+    v2c = v0d.copy()
+    v2a.rotate(q2)
+    v2b.rotate(q2)
+    v2c.rotate(q2)
+    v2a += pp2
+    v2b += pp2
+    v2c += pp2
+
+    pp3 = p0 + d3 / 3
+    v3a = v0a.copy()
+    v3b = v0b.copy()
+    v3c = v0c.copy()
+    v3a.rotate(q3)
+    v3b.rotate(q3)
+    v3c.rotate(q3)
+    v3a += pp3
+    v3b += pp3
+    v3c += pp3
+
+    v0a += p0
+    v0b += p0
+    v0c += p0
+    v0d += p0
+
+    verts = (v1a, v1b, v1c, v2a, v2b, v2c, v3a, v3b, v3c, v0a, v0b, v0c, v0d)
+    faces = (
+        #(1, 2, 12, 11),
+        #(9, 12, 2, 0),
+        #(11, 9, 0, 1),
+
+        #(5, 4, 10, 12),
+        #(4, 3, 9, 10),
+        #(3, 5, 12, 9),
+
+        (8, 7, 11, 10),
+        (7, 5, 9, 11),
+        (6, 8, 10, 9),
+
+        (10, 11, 12))
+
+    return verts, faces
+
 
 def bridgequads(aquad, bquad, verts):
-  "return faces, aloop, bloop"
-  ai,bi,_= min([(ai,bi,(verts[a]-verts[b]).length_squared) for ai,a in enumerate(aquad) for bi,b in enumerate(bquad)],key=lambda x:x[2])
-  n=len(aquad)
-  #print([(aquad[(ai+i)%n], aquad[(ai+i+1)%n], bquad[(bi+i+1)%n], bquad[(bi+i)%n]) for i in range(n)],"\n", [aquad[(ai+i)%n] for i in range(n)],"\n",  [aquad[(bi+i)%n] for i in range(n)])
-  
-  #print('bridgequads',aquad,bquad,ai,bi)
-  
-  return ([(aquad[(ai+i)%n], aquad[(ai+i+1)%n], bquad[(bi+i+1)%n], bquad[(bi+i)%n]) for i in range(n)], [aquad[(ai+i)%n] for i in range(n)],  [bquad[(bi+i)%n] for i in range(n)])
-
-def quadfork(p0,p1,p2,p3,r0,r1,r2,r3):
-  d1=p1-p0
-  d2=p2-p0
-  d3=p3-p0
-  a=(d3-d2).normalized()
-  n=d2.cross(d3).normalized()
-  pp1=p0+d1/3
-  pp2=p0+d2/3
-  pp3=p0+d3/3
-
-  v2a=pp2+( n+a)*r2
-  v2b=pp2+( n-a)*r2
-  v2c=pp2+(-n-a)*r2
-  v2d=pp2+(-n+a)*r2
-
-  v3a=pp3+( n+a)*r3
-  v3b=pp3+( n-a)*r3
-  v3c=pp3+(-n-a)*r3
-  v3d=pp3+(-n+a)*r3
-
-  a=d1.cross(n).normalized()
-  n=a.cross(d1).normalized()
-  v1a=pp1+( n+a)*r1
-  v1b=pp1+( n-a)*r1
-  v1c=pp1+(-n-a)*r1
-  v1d=pp1+(-n+a)*r1
-
-  #the top of the connecting block consist of two quads
-  v0a=p0+( n+a)*r0
-  v0b=p0+( n-a)*r0
-  v0c=p0+(-n-a)*r0
-  v0d=p0+(-n+a)*r0
-  v0ab=p0+n*r0
-  v0cd=p0-n*r0
-  #the bottom is a single quad (which means the front and back are 5gons)
-  d=d1.normalized()*r0*0.1
-  vb0a=v0a+d
-  vb0b=v0b+d
-  vb0c=v0c+d
-  vb0d=v0d+d
-
-  verts=[v1a,v1b,v1c,v1d,            # 0 1 2 3
-         v2a,v2b,v2c,v2d,            # 4 5 6 7
-         v3a,v3b,v3c,v3d,            # 8 9 10 11
-         v0a,v0ab,v0b,v0c,v0cd,v0d,  # 12 13 14 15 16 17
-         vb0a,vb0b,vb0c,vb0d]        # 18 19 20 21
-
-  faces=[(0,1,19,18),       # p1->p0 bottom
-         (1,2,20,19),
-         (2,3,21,20),
-         (3,0,18,21),
-
-                #(4,5,14,13),       # p2 -> p0 top right
-         #(5,6,15,14),
-         #(6,7,16,15),
-         #(7,4,13,16),
-                
-                (13,14,5,4),
-                (14,15,6,5),
-                (15,16,7,6),
-                (16,13,4,7),
-
-         #(8,9,13,12),       # p3 -> p0 top left
-         #(9,10,16,13),
-         #(10,11,17,16),
-         #(11,8,12,17),
-
-                (12,13,9,8),
-                (13,16,10,9),
-                (16,17,11,10),
-                (17,12,8,11),
-                
-         #(12,17,21,18),     # connecting block
-         #(14,15,20,19),
-         #(12,13,14,19,18),
-         #(15,16,17,21,20)]
-
-                (12,17,21,18),     # connecting block
-         (19,20,15,14),
-         (18,19,14,13,12),
-         (20,21,17,16,15)]
-
-  return verts, faces
\ No newline at end of file
+    "return faces,  aloop,  bloop"
+    ai, bi, _ = min([(ai, bi, (verts[a] - verts[b]).length_squared) for ai, a in enumerate(aquad) for bi, b in enumerate(bquad)], key=lambda x: x[2])
+    n = len(aquad)
+    #print([(aquad[(ai+i)%n],  aquad[(ai+i+1)%n],  bquad[(bi+i+1)%n],  bquad[(bi+i)%n]) for i in range(n)], "\n",  [aquad[(ai+i)%n] for i in range(n)], "\n",   [aquad[(bi+i)%n] for i in range(n)])
+
+    #print('bridgequads', aquad, bquad, ai, bi)
+
+    return ([(aquad[(ai + i) % n], aquad[(ai + i + 1) % n], bquad[(bi + i + 1) % n], bquad[(bi + i) % n]) for i in range(n)], [aquad[(ai + i) % n] for i in range(n)], [bquad[(bi + i) % n] for i in range(n)])
+
+
+def quadfork(p0, p1, p2, p3, r0, r1, r2, r3):
+    d1 = p1 - p0
+    d2 = p2 - p0
+    d3 = p3 - p0
+    a = (d3 - d2).normalized()
+    n = d2.cross(d3).normalized()
+    pp1 = p0 + d1 / 3
+    pp2 = p0 + d2 / 3
+    pp3 = p0 + d3 / 3
+
+    v2a = pp2 + (n + a) * r2
+    v2b = pp2 + (n - a) * r2
+    v2c = pp2 + (-n - a) * r2
+    v2d = pp2 + (-n + a) * r2
+
+    v3a = pp3 + (n + a) * r3
+    v3b = pp3 + (n - a) * r3
+    v3c = pp3 + (-n - a) * r3
+    v3d = pp3 + (-n + a) * r3
+
+    a = d1.cross(n).normalized()
+    n = a.cross(d1).normalized()
+    v1a = pp1 + (n + a) * r1
+    v1b = pp1 + (n - a) * r1
+    v1c = pp1 + (-n - a) * r1
+    v1d = pp1 + (-n + a) * r1
+
+    #the top of the connecting block consist of two quads
+    v0a = p0 + (n + a) * r0
+    v0b = p0 + (n - a) * r0
+    v0c = p0 + (-n - a) * r0
+    v0d = p0 + (-n + a) * r0
+    v0ab = p0 + n * r0
+    v0cd = p0 - n * r0
+    #the bottom is a single quad (which means the front and back are 5gons)
+    d = d1.normalized() * r0 * 0.1
+    vb0a = v0a + d
+    vb0b = v0b + d
+    vb0c = v0c + d
+    vb0d = v0d + d
+
+    verts = [v1a, v1b, v1c, v1d,             # 0 1 2 3
+        v2a, v2b, v2c, v2d,             # 4 5 6 7
+        v3a, v3b, v3c, v3d,             # 8 9 10 11
+        v0a, v0ab, v0b, v0c, v0cd, v0d,   # 12 13 14 15 16 17
+        vb0a, vb0b, vb0c, vb0d]        # 18 19 20 21
+
+    faces = [(0, 1, 19, 18),        # p1->p0 bottom
+        (1, 2, 20, 19),
+        (2, 3, 21, 20),
+        (3, 0, 18, 21),
+
+        #(4, 5, 14, 13),        # p2 -> p0 top right
+        #(5, 6, 15, 14),
+        #(6, 7, 16, 15),
+        #(7, 4, 13, 16),
+
+        (13, 14, 5, 4),
+        (14, 15, 6, 5),
+        (15, 16, 7, 6),
+        (16, 13, 4, 7),
+
+        #(8, 9, 13, 12),        # p3 -> p0 top left
+        #(9, 10, 16, 13),
+        #(10, 11, 17, 16),
+        #(11, 8, 12, 17),
+
+        (12, 13, 9, 8),
+        (13, 16, 10, 9),
+        (16, 17, 11, 10),
+        (17, 12, 8, 11),
+
+        #(12, 17, 21, 18),      # connecting block
+        #(14, 15, 20, 19),
+        #(12, 13, 14, 19, 18),
+        #(15, 16, 17, 21, 20)]
+
+        (12, 17, 21, 18),      # connecting block
+        (19, 20, 15, 14),
+        (18, 19, 14, 13, 12),
+        (20, 21, 17, 16, 15)]
+
+    return verts, faces
index 326e1fe..8aa18de 100644 (file)
@@ -1,22 +1,24 @@
+# <pep8 compliant>
 from collections import OrderedDict
 from time import time
 
+
 class Timer:
-       "an ordered collection of timestamps."
+    """an ordered collection of timestamps."""
 
-       def __init__(self):
-               self.od = OrderedDict()
-               self.od['Start']=time()
+    def __init__(self):
+        self.od = OrderedDict()
+        self.od['Start'] = time()
 
-       def add(self,label):
-               "add a new labeled timestamp"
-               self.od[label]=time()
+    def add(self, label):
+        """add a new labeled timestamp"""
+        self.od[label] = time()
 
-       def __str__(self):
-               "print a list of timings in order of addition. Second column is time since start."
-               keys=list(self.od.keys())
-               if len(keys)<2:
-                       return ""
-               return "\n".join("%-20s: %6.1fs %6.1f"%(keys[i],
-                       self.od[keys[i]]-self.od[keys[i-1]],
-                       self.od[keys[i]]-self.od[keys[0]]) for i in range(1,len(keys)))
+    def __str__(self):
+        """print a list of timings in order of addition. Second column is time since start."""
+        keys = list(self.od.keys())
+        if len(keys) < 2:
+            return ""
+        return "\n".join("%-20s: %6.1fs %6.1f" % (keys[i],
+            self.od[keys[i]] - self.od[keys[i - 1]],
+            self.od[keys[i]] - self.od[keys[0]]) for i in range(1, len(keys)))