Added a 2d convex hull function to BPyMathutils
[blender.git] / release / scripts / bpymodules / BPyMathutils.py
1 # $Id$
2 #
3 # --------------------------------------------------------------------------
4 # helper functions to be used by other scripts
5 # --------------------------------------------------------------------------
6 # ***** BEGIN GPL LICENSE BLOCK *****
7 #
8 # This program is free software; you can redistribute it and/or
9 # modify it under the terms of the GNU General Public License
10 # as published by the Free Software Foundation; either version 2
11 # of the License, or (at your option) any later version.
12 #
13 # This program is distributed in the hope that it will be useful,
14 # but WITHOUT ANY WARRANTY; without even the implied warranty of
15 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16 # GNU General Public License for more details.
17 #
18 # You should have received a copy of the GNU General Public License
19 # along with this program; if not, write to the Free Software Foundation,
20 # Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
21 #
22 # ***** END GPL LICENCE BLOCK *****
23 # --------------------------------------------------------------------------
24
25 import Blender
26 from Blender.Mathutils import *
27
28 # ------ Mersenne Twister - start
29
30 # Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura.
31 # Any feedback is very welcome. For any question, comments,
32 # see http://www.math.keio.ac.jp/matumoto/emt.html or email
33 # matumoto@math.keio.ac.jp
34
35 # The link above is dead, this is the new one:
36 # http://www.math.sci.hiroshima-u.ac.jp/m-mat/MT/emt.html
37 # And here the license info, from Mr. Matsumoto's site:
38 # Until 2001/4/6, MT had been distributed under GNU Public License,
39 # but after 2001/4/6, we decided to let MT be used for any purpose, including
40 # commercial use. 2002-versions mt19937ar.c, mt19937ar-cok.c are considered
41 # to be usable freely.
42 #
43 # So from the year above (1997), this code is under GPL.
44
45 # Period parameters
46 N = 624
47 M = 397
48 MATRIX_A = 0x9908b0dfL   # constant vector a
49 UPPER_MASK = 0x80000000L # most significant w-r bits
50 LOWER_MASK = 0x7fffffffL # least significant r bits
51
52 # Tempering parameters
53 TEMPERING_MASK_B = 0x9d2c5680L
54 TEMPERING_MASK_C = 0xefc60000L
55
56 def TEMPERING_SHIFT_U(y):
57     return (y >> 11)
58
59 def TEMPERING_SHIFT_S(y):
60     return (y << 7)
61
62 def TEMPERING_SHIFT_T(y):
63     return (y << 15)
64
65 def TEMPERING_SHIFT_L(y):
66     return (y >> 18)
67
68 mt = []   # the array for the state vector
69 mti = N+1 # mti==N+1 means mt[N] is not initialized
70
71 # initializing the array with a NONZERO seed
72 def sgenrand(seed):
73   # setting initial seeds to mt[N] using
74   # the generator Line 25 of Table 1 in
75   # [KNUTH 1981, The Art of Computer Programming
76   #    Vol. 2 (2nd Ed.), pp102]
77
78   global mt, mti
79
80   mt = []
81   
82   mt.append(seed & 0xffffffffL)
83   for i in xrange(1, N + 1):
84     mt.append((69069 * mt[i-1]) & 0xffffffffL)
85
86   mti = i
87 # end sgenrand
88
89
90 def genrand():
91   global mt, mti
92   
93   mag01 = [0x0L, MATRIX_A]
94   # mag01[x] = x * MATRIX_A  for x=0,1
95   y = 0
96   
97   if mti >= N: # generate N words at one time
98     if mti == N+1:   # if sgenrand() has not been called,
99       sgenrand(4357) # a default initial seed is used
100
101     for kk in xrange((N-M) + 1):
102       y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK)
103       mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1]
104
105     for kk in xrange(kk, N):
106       y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK)
107       mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1]
108
109     y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK)
110     mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1]
111
112     mti = 0
113
114   y = mt[mti]
115   mti += 1
116   y ^= TEMPERING_SHIFT_U(y)
117   y ^= TEMPERING_SHIFT_S(y) & TEMPERING_MASK_B
118   y ^= TEMPERING_SHIFT_T(y) & TEMPERING_MASK_C
119   y ^= TEMPERING_SHIFT_L(y)
120
121   return ( float(y) / 0xffffffffL ) # reals
122
123 #------ Mersenne Twister -- end 
124
125
126
127
128 """ 2d convexhull
129 Based from Dinu C. Gherman's work,
130 modified for Blender/Mathutils by Campell Barton
131 """
132 ######################################################################
133 # Public interface
134 ######################################################################
135 from Blender.Mathutils import DotVecs
136 def convexHull(point_list_2d):
137         """Calculate the convex hull of a set of vectors
138         The vectors can be 3 or 4d but only the Xand Y are used.
139         returns a list of convex hull indicies to the given point list
140         """
141
142         ######################################################################
143         # Helpers
144         ######################################################################
145
146         def _myDet(p, q, r):
147                 """Calc. determinant of a special matrix with three 2D points.
148
149                 The sign, "-" or "+", determines the side, right or left,
150                 respectivly, on which the point r lies, when measured against
151                 a directed vector from p to q.
152                 """
153                 return (q.x*r.y + p.x*q.y + r.x*p.y)  -  (q.x*p.y + r.x*q.y + p.x*r.y)
154         
155         def _isRightTurn((p, q, r)):
156                 "Do the vectors pq:qr form a right turn, or not?"
157                 #assert p[0] != q[0] and q[0] != r[0] and p[0] != r[0]
158                 if _myDet(p[0], q[0], r[0]) < 0:
159                         return 1
160                 else:
161                         return 0
162
163         # Get a local list copy of the points and sort them lexically.
164         points = [(p, i) for i, p in enumerate(point_list_2d)]
165         points.sort(lambda a,b: cmp((a[0].x, a[0].y), (b[0].x, b[0].y)))
166
167         # Build upper half of the hull.
168         upper = [points[0], points[1]] # cant remove these.
169         for i in xrange(len(points)-2):
170                 upper.append(points[i+2])
171                 while len(upper) > 2 and not _isRightTurn(upper[-3:]):
172                         del upper[-2]
173
174         # Build lower half of the hull.
175         points.reverse()
176         lower = [points.pop(0), points.pop(1)]
177         for p in points:
178                 lower.append(p)
179                 while len(lower) > 2 and not _isRightTurn(lower[-3:]):
180                         del lower[-2]
181         
182         # Concatenate both halfs and return.
183         return [p[1] for ls in (upper, lower) for p in ls]
184
185
186 def lineIntersect2D(v1a, v1b,  v2a, v2b):
187         '''
188         Do 2 lines intersect, if so where.
189         If there is an error, the retured X value will be None
190         the y will be an error code- usefull when debugging.
191         
192         the first line is (v1a, v1b)
193         the second is (v2a, v2b)
194         by Campbell Barton
195         This function accounts for all known cases of 2 lines ;)
196         '''
197         
198         x1,y1=          v1a.x, v1a.y
199         x2,y2=          v1b.x, v1b.y
200         _x1,_y1=        v2a.x, v2a.y
201         _x2,_y2=        v2b.x, v2b.y
202
203         # Bounding box intersection first.
204         if min(x1, x2) > max(_x1, _x2) or \
205         max(x1, x2) < min(_x1, _x2) or \
206         min(y1, y2) > max(_y1, _y2) or \
207         max(y1, y2) < min(_y1, _y2):
208                 return None, 100 # Basic Bounds intersection TEST returns false.
209         
210         # are either of the segments points? Check Seg1
211         if abs(x1 - x2) + abs(y1 - y2) <= SMALL_NUM:
212                 return None, 101
213         
214         # are either of the segments points? Check Seg2
215         if abs(_x1 - _x2) + abs(_y1 - _y2) <= SMALL_NUM:
216                 return None, 102
217         
218         # Make sure the HOZ/Vert Line Comes first.
219         if abs(_x1 - _x2) < SMALL_NUM or abs(_y1 - _y2) < SMALL_NUM:
220                 x1, x2, y1, y2, _x1, _x2, _y1, _y2 = _x1, _x2, _y1, _y2, x1, x2, y1, y2
221         
222         if abs(x2-x1) < SMALL_NUM: # VERTICLE LINE
223                 if abs(_x2-_x1) < SMALL_NUM: # VERTICLE LINE SEG2
224                         return None, 111 # 2 verticle lines dont intersect.
225                 
226                 elif abs(_y2-_y1) < SMALL_NUM:
227                         return x1, _y1 # X of vert, Y of hoz. no calculation.           
228                 
229                 yi = ((_y1 / abs(_x1 - _x2)) * abs(_x2 - x1)) + ((_y2 / abs(_x1 - _x2)) * abs(_x1 - x1))
230                 
231                 if yi > max(y1, y2): # New point above seg1's vert line
232                         return None, 112
233                 elif yi < min(y1, y2): # New point below seg1's vert line
234                         return None, 113
235                         
236                 return x1, yi # Intersecting.
237         
238         
239         if abs(y2-y1) < SMALL_NUM: # HOZ LINE
240                 if abs(_y2-_y1) < SMALL_NUM: # HOZ LINE SEG2
241                         return None, 121 # 2 hoz lines dont intersect.
242                 
243                 # Can skip vert line check for seg 2 since its covered above.   
244                 xi = ((_x1 / abs(_y1 - _y2)) * abs(_y2 - y1)) + ((_x2 / abs(_y1 - _y2)) * abs(_y1 - y1))
245                 if xi > max(x1, x2): # New point right of seg1's hoz line
246                         return None, 112
247                 elif xi < min(x1, x2): # New point left of seg1's hoz line
248                         return None, 113
249                 
250                 return xi, y1 # Intersecting.
251         
252         # Accounted for hoz/vert lines. Go on with both anglular.
253         b1 = (y2-y1)/(x2-x1)
254         b2 = (_y2-_y1)/(_x2-_x1)
255         a1 = y1-b1*x1
256         a2 = _y1-b2*_x1
257         
258         if b1 - b2 == 0.0:
259                 return None, 200        
260         
261         xi = - (a1-a2)/(b1-b2)
262         yi = a1+b1*xi
263         if (x1-xi)*(xi-x2) >= 0 and (_x1-xi)*(xi-_x2) >= 0 and (y1-yi)*(yi-y2) >= 0 and (_y1-yi)*(yi-_y2)>=0:
264                 return xi, yi
265         else:
266                 return None, 300
267
268