7ebaefa3017c27a18129cd6ecb8987581c1eb88d
[synfig.git] / synfig-core / trunk / src / modules / mod_noise / random.cpp
1 /* === S Y N F I G ========================================================= */
2 /*!     \file mod_noise/random.cpp
3 **      \brief blehh
4 **
5 **      $Id$
6 **
7 **      \legal
8 **      Copyright (c) 2002-2005 Robert B. Quattlebaum Jr., Adrian Bentley
9 **
10 **      This package is free software; you can redistribute it and/or
11 **      modify it under the terms of the GNU General Public License as
12 **      published by the Free Software Foundation; either version 2 of
13 **      the License, or (at your option) any later version.
14 **
15 **      This package is distributed in the hope that it will be useful,
16 **      but WITHOUT ANY WARRANTY; without even the implied warranty of
17 **      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 **      General Public License for more details.
19 **      \endlegal
20 */
21 /* ========================================================================= */
22
23 /* === H E A D E R S ======================================================= */
24
25 #ifdef USING_PCH
26 #       include "pch.h"
27 #else
28 #ifdef HAVE_CONFIG_H
29 #       include <config.h>
30 #endif
31
32 #include "random.h"
33 #include <cmath>
34 #include <cstdlib>
35 #endif
36
37 // A fast 32-bit linear congruential random number generator
38 class quick_rng
39 {
40         unsigned long next;
41 public:
42         quick_rng(unsigned long seed):next(seed) { }
43
44         void set_seed(unsigned long x)
45         {
46                 next=x;
47         }
48
49         unsigned long i32()
50         {
51                 static const unsigned long a(1664525);
52                 static const unsigned long c(1013904223);
53
54                 return next=next*a+c;
55         }
56
57         unsigned long i16()
58         {
59                 return i32()>>16;
60         }
61
62         float f()
63         {
64                 static const float m(int(65535));
65
66                 return float(i16())/m;
67         }
68 };
69
70 /* === M A C R O S ========================================================= */
71
72 #define PI      (3.1415927)
73
74 /* === G L O B A L S ======================================================= */
75
76 /* === P R O C E D U R E S ================================================= */
77
78 /* === M E T H O D S ======================================================= */
79
80 void
81 Random::set_seed(int x)
82 {
83         seed_=x;
84 }
85
86 float
87 Random::operator()(const int salt,const int x,const int y,const int t)const
88 {
89         static const unsigned int a(21870);
90         static const unsigned int b(11213);
91         static const unsigned int c(36979);
92         static const unsigned int d(31337);
93
94         quick_rng rng(
95                 ( static_cast<unsigned int>(x+y)        * a ) ^
96                 ( static_cast<unsigned int>(y+t)        * b ) ^
97                 ( static_cast<unsigned int>(t+x)        * c ) ^
98                 ( static_cast<unsigned int>(seed_+salt) * d )
99         );
100
101         return rng.f() * 2.0f - 1.0f;
102 }
103
104 float
105 Random::operator()(int smooth,int subseed,float xf,float yf,float tf)const
106 {
107         int x((int)floor(xf));
108         int y((int)floor(yf));
109         int t((int)floor(tf));
110
111         switch(smooth)
112         {
113         case 4: // cubic
114                 {
115                         #define f(j,i,k)        ((*this)(subseed,i,j,k))
116                         //Using catmull rom interpolation because it doesn't blur at all
117                         //bezier curve with intermediate ctrl pts: 0.5/3(p(i+1) - p(i-1)) and similar
118                         float xfa [4], tfa[4];
119
120                         //precalculate indices (all clamped) and offset
121                         const int xa[] = {x-1,x,x+1,x+2};
122
123                         const int ya[] = {y-1,y,y+1,y+2};
124
125                         const int ta[] = {t-1,t,t+1,t+2};
126
127                         const float dx(xf-x);
128                         const float dy(yf-y);
129                         const float dt(tf-t);
130
131                         //figure polynomials for each point
132                         const float txf[] =
133                         {
134                                 0.5*dx*(dx*(dx*(-1) + 2) - 1),  //-t + 2t^2 -t^3
135                                 0.5*(dx*(dx*(3*dx - 5)) + 2),   //2 - 5t^2 + 3t^3
136                                 0.5*dx*(dx*(-3*dx + 4) + 1),    //t + 4t^2 - 3t^3
137                                 0.5*dx*dx*(dx-1)                                //-t^2 + t^3
138                         };
139
140                         const float tyf[] =
141                         {
142                                 0.5*dy*(dy*(dy*(-1) + 2) - 1),  //-t + 2t^2 -t^3
143                                 0.5*(dy*(dy*(3*dy - 5)) + 2),   //2 - 5t^2 + 3t^3
144                                 0.5*dy*(dy*(-3*dy + 4) + 1),    //t + 4t^2 - 3t^3
145                                 0.5*dy*dy*(dy-1)                                //-t^2 + t^3
146                         };
147
148                         const float ttf[] =
149                         {
150                                 0.5*dt*(dt*(dt*(-1) + 2) - 1),  //-t + 2t^2 -t^3
151                                 0.5*(dt*(dt*(3*dt - 5)) + 2),   //2 - 5t^2 + 3t^3
152                                 0.5*dt*(dt*(-3*dt + 4) + 1),    //t + 4t^2 - 3t^3
153                                 0.5*dt*dt*(dt-1)                                //-t^2 + t^3
154                         };
155
156                         //evaluate polynomial for each row
157                         for(int i = 0; i < 4; ++i)
158                         {
159                                 for(int j = 0; j < 4; ++j)
160                                 {
161                                         tfa[j] = f(ya[i],xa[j],ta[0])*ttf[0] + f(ya[i],xa[j],ta[1])*ttf[1] + f(ya[i],xa[j],ta[2])*ttf[2] + f(ya[i],xa[j],ta[3])*ttf[3];
162                                 }
163                                 xfa[i] = tfa[0]*txf[0] + tfa[1]*txf[1] + tfa[2]*txf[2] + tfa[3]*txf[3];
164                         }
165
166                         //return the cumulative column evaluation
167                         return xfa[0]*tyf[0] + xfa[1]*tyf[1] + xfa[2]*tyf[2] + xfa[3]*tyf[3];
168 #undef f
169                 }
170                 break;
171
172
173         case 5: // Fast Spline (non-animated)
174                 {
175 #define P(x)    (((x)>0)?((x)*(x)*(x)):0.0f)
176 #define R(x)    ( P(x+2) - 4.0f*P(x+1) + 6.0f*P(x) - 4.0f*P(x-1) )*(1.0f/6.0f)
177 #define F(i,j)  ((*this)(subseed,i+x,j+y)*(R((i)-a)*R(b-(j))))
178 #define FT(i,j,k)       ((*this)(subseed,i+x,j+y,k+t)*(R((i)-a)*R(b-(j))*R((k)-c)))
179 #define Z(i,j) ret+=F(i,j)
180 #define ZT(i,j,k) ret+=FT(i,j,k)
181 #define X(i,j)  // placeholder... To make box more symmetric
182 #define XT(i,j,k)       // placeholder... To make box more symmetric
183
184                 float a(xf-x), b(yf-y);
185
186                 // Interpolate
187                 float ret(F(0,0));
188                 Z(-1,-1); Z(-1, 0); Z(-1, 1); Z(-1, 2);
189                 Z( 0,-1); X( 0, 0); Z( 0, 1); Z( 0, 2);
190                 Z( 1,-1); Z( 1, 0); Z( 1, 1); Z( 1, 2);
191                 Z( 2,-1); Z( 2, 0); Z( 2, 1); Z( 2, 2);
192
193                 return ret;
194         }
195
196         case 3: // Spline (animated)
197                 {
198                         float a(xf-x), b(yf-y), c(tf-t);
199
200                         // Interpolate
201                         float ret(FT(0,0,0));
202                         ZT(-1,-1,-1); ZT(-1, 0,-1); ZT(-1, 1,-1); ZT(-1, 2,-1);
203                         ZT( 0,-1,-1); ZT( 0, 0,-1); ZT( 0, 1,-1); ZT( 0, 2,-1);
204                         ZT( 1,-1,-1); ZT( 1, 0,-1); ZT( 1, 1,-1); ZT( 1, 2,-1);
205                         ZT( 2,-1,-1); ZT( 2, 0,-1); ZT( 2, 1,-1); ZT( 2, 2,-1);
206
207                         ZT(-1,-1, 0); ZT(-1, 0, 0); ZT(-1, 1, 0); ZT(-1, 2, 0);
208                         ZT( 0,-1, 0); XT( 0, 0, 0); ZT( 0, 1, 0); ZT( 0, 2, 0);
209                         ZT( 1,-1, 0); ZT( 1, 0, 0); ZT( 1, 1, 0); ZT( 1, 2, 0);
210                         ZT( 2,-1, 0); ZT( 2, 0, 0); ZT( 2, 1, 0); ZT( 2, 2, 0);
211
212                         ZT(-1,-1, 1); ZT(-1, 0, 1); ZT(-1, 1, 1); ZT(-1, 2, 1);
213                         ZT( 0,-1, 1); ZT( 0, 0, 1); ZT( 0, 1, 1); ZT( 0, 2, 1);
214                         ZT( 1,-1, 1); ZT( 1, 0, 1); ZT( 1, 1, 1); ZT( 1, 2, 1);
215                         ZT( 2,-1, 1); ZT( 2, 0, 1); ZT( 2, 1, 1); ZT( 2, 2, 1);
216
217                         ZT(-1,-1, 2); ZT(-1, 0, 2); ZT(-1, 1, 2); ZT(-1, 2, 2);
218                         ZT( 0,-1, 2); ZT( 0, 0, 2); ZT( 0, 1, 2); ZT( 0, 2, 2);
219                         ZT( 1,-1, 2); ZT( 1, 0, 2); ZT( 1, 1, 2); ZT( 1, 2, 2);
220                         ZT( 2,-1, 2); ZT( 2, 0, 2); ZT( 2, 1, 2); ZT( 2, 2, 2);
221
222                         return ret;
223
224 /*
225
226                         float dx=xf-x;
227                         float dy=yf-y;
228                         float dt=tf-t;
229
230                         float ret=0;
231                         int i,j,h;
232                         for(h=-1;h<=2;h++)
233                                 for(i=-1;i<=2;i++)
234                                         for(j=-1;j<=2;j++)
235                                                 ret+=(*this)(subseed,i+x,j+y,h+t)*(R(i-dx)*R(j-dy)*R(h-dt));
236                         return ret;
237 */
238                 }
239                 break;
240 #undef X
241 #undef Z
242 #undef F
243 #undef P
244 #undef R
245
246         case 2: // Cosine
247         if((float)t==tf)
248         {
249                 int x((int)floor(xf));
250                 int y((int)floor(yf));
251                 float a=xf-x;
252                 float b=yf-y;
253                 a=(1.0f-cos(a*PI))*0.5f;
254                 b=(1.0f-cos(b*PI))*0.5f;
255                 float c=1.0-a;
256                 float d=1.0-b;
257                 int x2=x+1,y2=y+1;
258                 return
259                         (*this)(subseed,x,y,t)*(c*d)+
260                         (*this)(subseed,x2,y,t)*(a*d)+
261                         (*this)(subseed,x,y2,t)*(c*b)+
262                         (*this)(subseed,x2,y2,t)*(a*b);
263         }
264         else
265         {
266                 float a=xf-x;
267                 float b=yf-y;
268                 float c=tf-t;
269
270                 a=(1.0f-cos(a*3.1415927))*0.5f;
271                 b=(1.0f-cos(b*3.1415927))*0.5f;
272
273                 // We don't perform this on the time axis, otherwise we won't
274                 // get smooth motion
275                 //c=(1.0f-cos(c*3.1415927))*0.5f;
276
277                 float d=1.0-a;
278                 float e=1.0-b;
279                 float f=1.0-c;
280
281                 int x2=x+1,y2=y+1,t2=t+1;
282
283                 return
284                         (*this)(subseed,x,y,t)*(d*e*f)+
285                         (*this)(subseed,x2,y,t)*(a*e*f)+
286                         (*this)(subseed,x,y2,t)*(d*b*f)+
287                         (*this)(subseed,x2,y2,t)*(a*b*f)+
288                         (*this)(subseed,x,y,t2)*(d*e*c)+
289                         (*this)(subseed,x2,y,t2)*(a*e*c)+
290                         (*this)(subseed,x,y2,t2)*(d*b*c)+
291                         (*this)(subseed,x2,y2,t2)*(a*b*c);
292         }
293         case 1: // Linear
294         if((float)t==tf)
295         {
296                 int x((int)floor(xf));
297                 int y((int)floor(yf));
298                 float a=xf-x;
299                 float b=yf-y;
300                 float c=1.0-a;
301                 float d=1.0-b;
302                 int x2=x+1,y2=y+1;
303                 return
304                         (*this)(subseed,x,y,t)*(c*d)+
305                         (*this)(subseed,x2,y,t)*(a*d)+
306                         (*this)(subseed,x,y2,t)*(c*b)+
307                         (*this)(subseed,x2,y2,t)*(a*b);
308         }
309         else
310         {
311
312                 float a=xf-x;
313                 float b=yf-y;
314                 float c=tf-t;
315
316                 float d=1.0-a;
317                 float e=1.0-b;
318                 float f=1.0-c;
319
320                 int x2=x+1,y2=y+1,t2=t+1;
321
322                 return
323                         (*this)(subseed,x,y,t)*(d*e*f)+
324                         (*this)(subseed,x2,y,t)*(a*e*f)+
325                         (*this)(subseed,x,y2,t)*(d*b*f)+
326                         (*this)(subseed,x2,y2,t)*(a*b*f)+
327                         (*this)(subseed,x,y,t2)*(d*e*c)+
328                         (*this)(subseed,x2,y,t2)*(a*e*c)+
329                         (*this)(subseed,x,y2,t2)*(d*b*c)+
330                         (*this)(subseed,x2,y2,t2)*(a*b*c);
331         }
332         default:
333         case 0:
334                 return (*this)(subseed,x,y,t);
335         }
336 }