Integrated Freestyle to rendering pipeline
[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 def convexHull(point_list_2d):
136         """Calculate the convex hull of a set of vectors
137         The vectors can be 3 or 4d but only the Xand Y are used.
138         returns a list of convex hull indicies to the given point list
139         """
140
141         ######################################################################
142         # Helpers
143         ######################################################################
144
145         def _myDet(p, q, r):
146                 """Calc. determinant of a special matrix with three 2D points.
147
148                 The sign, "-" or "+", determines the side, right or left,
149                 respectivly, on which the point r lies, when measured against
150                 a directed vector from p to q.
151                 """
152                 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)
153         
154         def _isRightTurn((p, q, r)):
155                 "Do the vectors pq:qr form a right turn, or not?"
156                 #assert p[0] != q[0] and q[0] != r[0] and p[0] != r[0]
157                 if _myDet(p[0], q[0], r[0]) < 0:
158                         return 1
159                 else:
160                         return 0
161
162         # Get a local list copy of the points and sort them lexically.
163         points = [(p, i) for i, p in enumerate(point_list_2d)]
164         
165         try:    points.sort(key = lambda a: (a[0].x, a[0].y))
166         except: points.sort(lambda a,b: cmp((a[0].x, a[0].y), (b[0].x, b[0].y)))
167
168         # Build upper half of the hull.
169         upper = [points[0], points[1]] # cant remove these.
170         for i in xrange(len(points)-2):
171                 upper.append(points[i+2])
172                 while len(upper) > 2 and not _isRightTurn(upper[-3:]):
173                         del upper[-2]
174
175         # Build lower half of the hull.
176         points.reverse()
177         lower = [points.pop(0), points.pop(1)]
178         for p in points:
179                 lower.append(p)
180                 while len(lower) > 2 and not _isRightTurn(lower[-3:]):
181                         del lower[-2]
182         
183         # Concatenate both halfs and return.
184         return [p[1] for ls in (upper, lower) for p in ls]
185
186
187 def plane2mat(plane, normalize= False):
188         '''
189         Takes a plane and converts to a matrix
190         points between 0 and 1 are up
191         1 and 2 are right
192         assumes the plane has 90d corners
193         '''
194         cent= (plane[0]+plane[1]+plane[2]+plane[3] ) /4.0
195
196         
197         up= cent - ((plane[0]+plane[1])/2.0)
198         right= cent - ((plane[1]+plane[2])/2.0)
199         z= up.cross(right)
200         
201         if normalize:
202                 up.normalize()
203                 right.normalize()
204                 z.normalize()
205         
206         mat= Matrix(up, right, z)
207         
208         # translate
209         mat.resize4x4()
210         tmat= Blender.Mathutils.TranslationMatrix(cent)
211         return mat * tmat
212
213
214 # Used for mesh_solidify.py and mesh_wire.py
215
216 # returns a length from an angle
217 # Imaging a 2d space.
218 # there is a hoz line at Y1 going to inf on both X ends, never moves (LINEA)
219 # down at Y0 is a unit length line point up at (angle) from X0,Y0 (LINEB)
220 # This function returns the length of LINEB at the point it would intersect LINEA
221 # - Use this for working out how long to make the vector - differencing it from surrounding faces,
222 # import math
223 from math import pi, sin, cos, sqrt
224
225 def angleToLength(angle):
226         # Alredy accounted for
227         if angle < 0.000001:    return 1.0
228         else:                                   return abs(1.0 / cos(pi*angle/180));