Fix T54414: Address changes to object.ray_cast, cleanup
[blender-addons-contrib.git] / add_mesh_space_tree / sca.py
1 # ##### BEGIN GPL LICENSE BLOCK #####
2 #
3 #  SCA Tree Generator, a Blender add-on
4 #  (c) 2013 Michel J. Anders (varkenvarken)
5 #
6 #  This program is free software; you can redistribute it and/or
7 #  modify it under the terms of the GNU General Public License
8 #  as published by the Free Software Foundation; either version 2
9 #  of the License, or (at your option) any later version.
10 #
11 #  This program is distributed in the hope that it will be useful,
12 #  but WITHOUT ANY WARRANTY; without even the implied warranty of
13 #  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14 #  GNU General Public License for more details.
15 #
16 #  You should have received a copy of the GNU General Public License
17 #  along with this program; if not, write to the Free Software Foundation,
18 #  Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
19 #
20 # ##### END GPL LICENSE BLOCK #####
21
22 # <pep8 compliant>
23
24 from collections import defaultdict as dd
25 from random import random, seed, expovariate
26 from math import sqrt
27 from functools import partial
28
29 from mathutils import Vector
30
31 from .kdtree import Tree
32
33
34 class Branchpoint:
35     def __init__(self, p, parent):
36         self.v = Vector(p)
37         self.parent = parent
38         self.connections = 0
39         self.apex = None
40         self.shoot = None
41         self.index = None
42
43
44 def sphere(r, p):
45     r2 = r * r
46     while True:
47         x = (random() * 2 - 1) * r
48         y = (random() * 2 - 1) * r
49         z = (random() * 2 - 1) * r
50         if x * x + y * y + z * z <= r2:
51             yield p + Vector((x, y, z))
52
53
54 class SCA:
55     def __init__(self, NENDPOINTS=100, d=0.3, NBP=2000, KILLDIST=5, INFLUENCE=15,
56                 SEED=42, volume=partial(sphere, 5, Vector((0, 0, 8))), TROPISM=0.0,
57                 exclude=lambda p: False,
58                 startingpoints=[]):
59         seed(SEED)
60         self.d = d
61         self.NBP = NBP
62         self.KILLDIST = KILLDIST * KILLDIST * d * d
63         self.INFLUENCE = INFLUENCE * INFLUENCE * d * d
64         self.TROPISM = TROPISM
65         self.exclude = exclude
66         self.endpoints = []
67
68         self.volumepoint = volume()
69         for i in range(NENDPOINTS):
70             self.endpoints.append(next(self.volumepoint))
71
72         self.branchpoints = [Branchpoint((0, 0, 0), None)] if len(startingpoints) == 0 else startingpoints
73
74     def iterate(self, newendpointsper1000=0, maxtime=0.0):  # maxtime still ignored for now
75
76         endpointsadded = 0.0
77         niterations = 0.0
78         newendpointsper1000 /= 1000.0
79         # time to the first new 'endpoint add event'
80         t = expovariate(newendpointsper1000) if newendpointsper1000 > 0.0 else 1
81
82         while self.NBP > 0 and (len(self.endpoints) > 0):
83             self.NBP -= 1
84             closestsendpoints = dd(list)
85
86             kill = set()
87
88             for ei, e in enumerate(self.endpoints):
89                 distance = None
90                 closestbp = None
91                 for bi, b in enumerate(self.branchpoints):
92                     ddd = b.v - e
93                     ddd = ddd.dot(ddd)
94                     if ddd < self.KILLDIST:
95                         kill.add(ei)
96                     elif (ddd < self.INFLUENCE and b.shoot is None) and ((distance is None) or (ddd < distance)):
97                         closestbp = bi
98                         distance = ddd
99                 if not (closestbp is None):
100                     closestsendpoints[closestbp].append(ei)
101
102             if len(closestsendpoints) < 1:
103                 break
104
105             for bi in closestsendpoints:
106                 sd = Vector((0, 0, 0))
107                 n = 0
108                 for ei in closestsendpoints[bi]:
109                     dv = self.branchpoints[bi].v - self.endpoints[ei]
110                     ll = sqrt(dv.dot(dv))
111                     sd -= dv / ll
112                     n += 1
113                 sd /= n
114                 ll = sqrt(sd.dot(sd))
115                 # don't know if this assumption is true:
116                 # if the unnormalised direction is very small, the endpoints are nearly
117                 # coplanar/colinear and at roughly the same distance
118                 # so no endpoints will be killed and we might end up adding the same branch again and again
119                 if ll < 1e-3:
120                     # print('SD very small')
121                     continue
122
123                 sd /= ll
124                 sd[2] += self.TROPISM
125                 ll = sqrt(sd.dot(sd))
126                 sd /= ll
127
128                 newp = self.branchpoints[bi].v + sd * self.d
129                 # the assumption we made earlier is not suffucient to prevent
130                 # adding the same branch so we need an expensive check:
131                 tooclose = False
132                 for dbi in self.branchpoints:
133                     dddd = newp - dbi.v
134                     if dddd.dot(dddd) < 1e-3:
135                         # print('BP to close to another')
136                         tooclose = True
137                 if tooclose:
138                     continue
139
140                 if not self.exclude(newp):
141                     bp = Branchpoint(newp, bi)
142                     self.branchpoints.append(bp)
143                     nbpi = len(self.branchpoints) - 1
144                     bp = self.branchpoints[bi]
145                     bp.connections += 1
146                     if bp.apex is None:
147                         bp.apex = nbpi
148                     else:
149                         bp.shoot = nbpi
150                     while not (bp.parent is None):
151                         bp = self.branchpoints[bp.parent]
152                         bp.connections += 1
153
154             self.endpoints = [ep for ei, ep in enumerate(self.endpoints) if not(ei in kill)]
155
156             if newendpointsper1000 > 0.0:
157                 # generate new endpoints with a poisson process
158                 # when we first arrive here, t already hold the time to the first event
159                 niterations += 1
160                 # we keep on adding endpoints as long as the next event still happens within this iteration
161                 while t < niterations:
162                     self.endpoints.append(next(self.volumepoint))
163                     endpointsadded += 1
164                     t += expovariate(newendpointsper1000)  # time to new 'endpoint add event'
165         """
166         if newendpointsper1000 > 0.0:
167             print("newendpoints / iteration %.3f, actual %.3f in %5.f iterations"%(
168                   newendpointsper1000, endpointsadded / niterations, niterations))
169         """
170     def iterate2(self, newendpointsper1000=0, maxtime=0.0):  # maxtime still ignored for now
171         """iterate using a kdtree fr the branchpoints"""
172         endpointsadded = 0.0
173         niterations = 0.0
174         newendpointsper1000 /= 1000.0
175         # time to the first new 'endpoint add event'
176         t = expovariate(newendpointsper1000) if newendpointsper1000 > 0.0 else 1
177
178         tree = Tree(3)
179         for bpi, bp in enumerate(self.branchpoints):
180             bp.index = bpi
181             tree.insert(bp.v, bp)
182
183         while self.NBP > 0 and (len(self.endpoints) > 0):
184             self.NBP -= 1
185             closestsendpoints = dd(list)
186
187             kill = set()
188
189             for ei, e in enumerate(self.endpoints):
190                 distance = None
191                 closestbp, distance = tree.nearest(e, checkempty=True)  # ignore forks, they have .data set to None
192                 if closestbp is not None:
193                     if distance < self.KILLDIST:
194                         kill.add(ei)
195                     elif distance < self.INFLUENCE:
196                         closestsendpoints[closestbp.data.index].append(ei)
197
198             if len(closestsendpoints) < 1:
199                 break
200
201             for bi in closestsendpoints:
202                 sd = Vector((0, 0, 0))
203                 n = 0
204                 for ei in closestsendpoints[bi]:
205                     dv = self.branchpoints[bi].v - self.endpoints[ei]
206                     dv.normalize()
207                     sd -= dv
208                     n += 1
209                 sd /= n
210                 ll = sd.length_squared
211                 # don't know if this assumption is true:
212                 # if the unnormalised direction is very small, the endpoints
213                 # are nearly coplanar/colinear and at roughly the same distance
214                 # so no endpoints will be killed and we might end up adding the same branch again and again
215                 if ll < 1e-3:
216                     # print('SD very small')
217                     continue
218
219                 sd /= ll
220                 sd[2] += self.TROPISM
221                 sd.normalize()
222
223                 newp = self.branchpoints[bi].v + sd * self.d
224                 # the assumption we made earlier is not suffucient to prevent
225                 # adding the same branch so we need an expensive check:
226                 _, dddd = tree.nearest(newp)  # no checkempty here, we want to check for forks as well
227                 if dddd < 1e-3:
228                         # print('BP to close to another')
229                         continue
230
231                 if not self.exclude(newp):
232                     bp = Branchpoint(newp, bi)
233                     self.branchpoints.append(bp)
234                     nbpi = len(self.branchpoints) - 1
235                     bp.index = nbpi
236                     tree.insert(bp.v, bp)
237                     bp = self.branchpoints[bi]
238                     bp.connections += 1
239                     if bp.apex is None:
240                         bp.apex = nbpi
241                     else:
242                         bp.shoot = nbpi
243                         node, _ = tree.nearest(bp.v)
244                         # signal that this is a fork so that we might ignore it
245                         # when searching for nearest neighbors
246                         node.data = None
247                     while not (bp.parent is None):
248                         bp = self.branchpoints[bp.parent]
249                         bp.connections += 1
250
251             self.endpoints = [ep for ei, ep in enumerate(self.endpoints) if not(ei in kill)]
252
253             if newendpointsper1000 > 0.0:
254                 # generate new endpoints with a poisson process
255                 # when we first arrive here, t already hold the time to the first event
256                 niterations += 1
257                 # we keep on adding endpoints as long as the next event still happens within this iteration
258                 while t < niterations:
259                     self.endpoints.append(next(self.volumepoint))
260                     endpointsadded += 1
261                     t += expovariate(newendpointsper1000)  # time to new 'endpoint add event'