synfig-core/trunk/src/modules/mod_particle/random.cpp

   1 /* === S Y N F I G ========================================================= */
   2 /*!     \file mod_particle/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
  36 #endif
  37
  38 /* === M A C R O S ========================================================= */
  39
  40 /* === G L O B A L S ======================================================= */
  41
  42 /* === P R O C E D U R E S ================================================= */
  43
  44 /* === M E T H O D S ======================================================= */
  45
  46 void
  47 Random::set_seed(int x)
  48 {
  49         seed_=x;
  50         srand(x);
  51         int i;
  52         for(i=0;i<POOL_SIZE;i++)
  53                 pool_[i]=rand();
  54
  55         x_mask=rand()+rand()*RAND_MAX;
  56         y_mask=rand()+rand()*RAND_MAX;
  57         t_mask=rand()+rand()*RAND_MAX;
  58 }
  59
  60 // this picks one of the POOL_SIZE (256) preset values out of the pool
  61 // and scales it to be in the range (-1, 1).  is that what it was
  62 // intended to do?  the distribution is pretty terrible, too, with
  63 // some elements being picked a hundred times more often than others
  64 float
  65 Random::operator()(const int salt,const int x,const int y,const int t)const
  66 {
  67         const int salt_hash(pool_[salt&(POOL_SIZE-1)]);
  68
  69         int index(((x^x_mask)+(y^y_mask)*234672+(t^t_mask)*8439573)^salt_hash);
  70
  71         index+=index*(index/POOL_SIZE);
  72
  73         return (float(pool_[index&(POOL_SIZE-1)])/float(RAND_MAX))*2.0f-1.0f;
  74 }
  75
  76 float
  77 Random::operator()(SmoothType smooth,int subseed,float xf,float yf,float tf)const
  78 {
  79         int x((int)floor(xf));
  80         int y((int)floor(yf));
  81         int t((int)floor(tf));
  82
  83         switch(smooth)
  84         {
  85         case SMOOTH_CUBIC:      // cubic
  86                 {
  87                         #define f(j,i,k)        ((*this)(subseed,i,j,k))
  88                         //Using catmull rom interpolation because it doesn't blur at all
  89                         // ( http://www.gamedev.net/reference/articles/article1497.asp )
  90                         //bezier curve with intermediate ctrl pts: 0.5/3(p(i+1) - p(i-1)) and similar
  91                         float xfa [4], tfa[4];
  92
  93                         //precalculate indices (all clamped) and offset
  94                         const int xa[] = {x-1,x,x+1,x+2};
  95
  96                         const int ya[] = {y-1,y,y+1,y+2};
  97
  98                         const int ta[] = {t-1,t,t+1,t+2};
  99
 100                         const float dx(xf-x);
 101                         const float dy(yf-y);
 102                         const float dt(tf-t);
 103
 104                         //figure polynomials for each point
 105                         const float txf[] =
 106                         {
 107                                 0.5*dx*(dx*(dx*(-1) + 2) - 1),  //-t + 2t^2 -t^3
 108                                 0.5*(dx*(dx*(3*dx - 5)) + 2),   //2 - 5t^2 + 3t^3
 109                                 0.5*dx*(dx*(-3*dx + 4) + 1),    //t + 4t^2 - 3t^3
 110                                 0.5*dx*dx*(dx-1)                                //-t^2 + t^3
 111                         };
 112
 113                         const float tyf[] =
 114                         {
 115                                 0.5*dy*(dy*(dy*(-1) + 2) - 1),  //-t + 2t^2 -t^3
 116                                 0.5*(dy*(dy*(3*dy - 5)) + 2),   //2 - 5t^2 + 3t^3
 117                                 0.5*dy*(dy*(-3*dy + 4) + 1),    //t + 4t^2 - 3t^3
 118                                 0.5*dy*dy*(dy-1)                                //-t^2 + t^3
 119                         };
 120
 121                         const float ttf[] =
 122                         {
 123                                 0.5*dt*(dt*(dt*(-1) + 2) - 1),  //-t + 2t^2 -t^3
 124                                 0.5*(dt*(dt*(3*dt - 5)) + 2),   //2 - 5t^2 + 3t^3
 125                                 0.5*dt*(dt*(-3*dt + 4) + 1),    //t + 4t^2 - 3t^3
 126                                 0.5*dt*dt*(dt-1)                                //-t^2 + t^3
 127                         };
 128
 129                         //evaluate polynomial for each row
 130                         for(int i = 0; i < 4; ++i)
 131                         {
 132                                 for(int j = 0; j < 4; ++j)
 133                                 {
 134                                         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];
 135                                 }
 136                                 xfa[i] = tfa[0]*txf[0] + tfa[1]*txf[1] + tfa[2]*txf[2] + tfa[3]*txf[3];
 137                         }
 138
 139                         //return the cumulative column evaluation
 140                         return xfa[0]*tyf[0] + xfa[1]*tyf[1] + xfa[2]*tyf[2] + xfa[3]*tyf[3];
 141 #undef f
 142                 }
 143                 break;
 144
 145
 146         case SMOOTH_FAST_SPLINE:        // Fast Spline (non-animated)
 147                 {
 148 #define P(x)    (((x)>0)?((x)*(x)*(x)):0.0f)
 149 #define R(x)    ( P(x+2) - 4.0f*P(x+1) + 6.0f*P(x) - 4.0f*P(x-1) )*(1.0f/6.0f)
 150 #define F(i,j)  ((*this)(subseed,i+x,j+y)*(R((i)-a)*R(b-(j))))
 151 #define FT(i,j,k)       ((*this)(subseed,i+x,j+y,k+t)*(R((i)-a)*R(b-(j))*R((k)-c)))
 152 #define Z(i,j) ret+=F(i,j)
 153 #define ZT(i,j,k) ret+=FT(i,j,k)
 154 #define X(i,j)  // placeholder... To make box more symmetric
 155 #define XT(i,j,k)       // placeholder... To make box more symmetric
 156
 157                 float a(xf-x), b(yf-y);
 158
 159                 // Interpolate
 160                 float ret(F(0,0));
 161                 Z(-1,-1); Z(-1, 0); Z(-1, 1); Z(-1, 2);
 162                 Z( 0,-1); X( 0, 0); Z( 0, 1); Z( 0, 2);
 163                 Z( 1,-1); Z( 1, 0); Z( 1, 1); Z( 1, 2);
 164                 Z( 2,-1); Z( 2, 0); Z( 2, 1); Z( 2, 2);
 165
 166                 return ret;
 167         }
 168
 169         case SMOOTH_SPLINE:     // Spline (animated)
 170                 {
 171                         float a(xf-x), b(yf-y), c(tf-t);
 172
 173                         // Interpolate
 174                         float ret(FT(0,0,0));
 175                         ZT(-1,-1,-1); ZT(-1, 0,-1); ZT(-1, 1,-1); ZT(-1, 2,-1);
 176                         ZT( 0,-1,-1); ZT( 0, 0,-1); ZT( 0, 1,-1); ZT( 0, 2,-1);
 177                         ZT( 1,-1,-1); ZT( 1, 0,-1); ZT( 1, 1,-1); ZT( 1, 2,-1);
 178                         ZT( 2,-1,-1); ZT( 2, 0,-1); ZT( 2, 1,-1); ZT( 2, 2,-1);
 179
 180                         ZT(-1,-1, 0); ZT(-1, 0, 0); ZT(-1, 1, 0); ZT(-1, 2, 0);
 181                         ZT( 0,-1, 0); XT( 0, 0, 0); ZT( 0, 1, 0); ZT( 0, 2, 0);
 182                         ZT( 1,-1, 0); ZT( 1, 0, 0); ZT( 1, 1, 0); ZT( 1, 2, 0);
 183                         ZT( 2,-1, 0); ZT( 2, 0, 0); ZT( 2, 1, 0); ZT( 2, 2, 0);
 184
 185                         ZT(-1,-1, 1); ZT(-1, 0, 1); ZT(-1, 1, 1); ZT(-1, 2, 1);
 186                         ZT( 0,-1, 1); ZT( 0, 0, 1); ZT( 0, 1, 1); ZT( 0, 2, 1);
 187                         ZT( 1,-1, 1); ZT( 1, 0, 1); ZT( 1, 1, 1); ZT( 1, 2, 1);
 188                         ZT( 2,-1, 1); ZT( 2, 0, 1); ZT( 2, 1, 1); ZT( 2, 2, 1);
 189
 190                         ZT(-1,-1, 2); ZT(-1, 0, 2); ZT(-1, 1, 2); ZT(-1, 2, 2);
 191                         ZT( 0,-1, 2); ZT( 0, 0, 2); ZT( 0, 1, 2); ZT( 0, 2, 2);
 192                         ZT( 1,-1, 2); ZT( 1, 0, 2); ZT( 1, 1, 2); ZT( 1, 2, 2);
 193                         ZT( 2,-1, 2); ZT( 2, 0, 2); ZT( 2, 1, 2); ZT( 2, 2, 2);
 194
 195                         return ret;
 196
 197 /*
 198
 199                         float dx=xf-x;
 200                         float dy=yf-y;
 201                         float dt=tf-t;
 202
 203                         float ret=0;
 204                         int i,j,h;
 205                         for(h=-1;h<=2;h++)
 206                                 for(i=-1;i<=2;i++)
 207                                         for(j=-1;j<=2;j++)
 208                                                 ret+=(*this)(subseed,i+x,j+y,h+t)*(R(i-dx)*R(j-dy)*R(h-dt));
 209                         return ret;
 210 */
 211                 }
 212                 break;
 213 #undef X
 214 #undef Z
 215 #undef F
 216 #undef P
 217 #undef R
 218
 219         case SMOOTH_COSINE:
 220         if((float)t==tf)
 221         {
 222                 int x((int)floor(xf));
 223                 int y((int)floor(yf));
 224                 float a=xf-x;
 225                 float b=yf-y;
 226                 a=(1.0f-cos(a*3.1415927))*0.5f;
 227                 b=(1.0f-cos(b*3.1415927))*0.5f;
 228                 float c=1.0-a;
 229                 float d=1.0-b;
 230                 int x2=x+1,y2=y+1;
 231                 return
 232                         (*this)(subseed,x,y,t)*(c*d)+
 233                         (*this)(subseed,x2,y,t)*(a*d)+
 234                         (*this)(subseed,x,y2,t)*(c*b)+
 235                         (*this)(subseed,x2,y2,t)*(a*b);
 236         }
 237         else
 238         {
 239                 float a=xf-x;
 240                 float b=yf-y;
 241                 float c=tf-t;
 242
 243                 a=(1.0f-cos(a*3.1415927))*0.5f;
 244                 b=(1.0f-cos(b*3.1415927))*0.5f;
 245
 246                 // We don't perform this on the time axis, otherwise we won't
 247                 // get smooth motion
 248                 //c=(1.0f-cos(c*3.1415927))*0.5f;
 249
 250                 float d=1.0-a;
 251                 float e=1.0-b;
 252                 float f=1.0-c;
 253
 254                 int x2=x+1,y2=y+1,t2=t+1;
 255
 256                 return
 257                         (*this)(subseed,x,y,t)*(d*e*f)+
 258                         (*this)(subseed,x2,y,t)*(a*e*f)+
 259                         (*this)(subseed,x,y2,t)*(d*b*f)+
 260                         (*this)(subseed,x2,y2,t)*(a*b*f)+
 261                         (*this)(subseed,x,y,t2)*(d*e*c)+
 262                         (*this)(subseed,x2,y,t2)*(a*e*c)+
 263                         (*this)(subseed,x,y2,t2)*(d*b*c)+
 264                         (*this)(subseed,x2,y2,t2)*(a*b*c);
 265         }
 266         case SMOOTH_LINEAR:
 267         if((float)t==tf)
 268         {
 269                 int x((int)floor(xf));
 270                 int y((int)floor(yf));
 271                 float a=xf-x;
 272                 float b=yf-y;
 273                 float c=1.0-a;
 274                 float d=1.0-b;
 275                 int x2=x+1,y2=y+1;
 276                 return
 277                         (*this)(subseed,x,y,t)*(c*d)+
 278                         (*this)(subseed,x2,y,t)*(a*d)+
 279                         (*this)(subseed,x,y2,t)*(c*b)+
 280                         (*this)(subseed,x2,y2,t)*(a*b);
 281         }
 282         else
 283         {
 284
 285                 float a=xf-x;
 286                 float b=yf-y;
 287                 float c=tf-t;
 288
 289                 float d=1.0-a;
 290                 float e=1.0-b;
 291                 float f=1.0-c;
 292
 293                 int x2=x+1,y2=y+1,t2=t+1;
 294
 295                 return
 296                         (*this)(subseed,x,y,t)*(d*e*f)+
 297                         (*this)(subseed,x2,y,t)*(a*e*f)+
 298                         (*this)(subseed,x,y2,t)*(d*b*f)+
 299                         (*this)(subseed,x2,y2,t)*(a*b*f)+
 300                         (*this)(subseed,x,y,t2)*(d*e*c)+
 301                         (*this)(subseed,x2,y,t2)*(a*e*c)+
 302                         (*this)(subseed,x,y2,t2)*(d*b*c)+
 303                         (*this)(subseed,x2,y2,t2)*(a*b*c);
 304         }
 305         default:
 306         case SMOOTH_DEFAULT:
 307                 return (*this)(subseed,x,y,t);
 308         }
 309 }