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

   1 /* === S I N F G =========================================================== */
   2 /*!     \file noise.cpp
   3 **      \brief blehh
   4 **
   5 **      $Id: random.cpp,v 1.1.1.1 2005/01/04 01:23:14 darco Exp $
   6 **
   7 **      \legal
   8 **      Copyright (c) 2002 Robert B. Quattlebaum Jr.
   9 **
  10 **      This software and associated documentation
  11 **      are CONFIDENTIAL and PROPRIETARY property of
  12 **      the above-mentioned copyright holder.
  13 **
  14 **      You may not copy, print, publish, or in any
  15 **      other way distribute this software without
  16 **      a prior written agreement with
  17 **      the copyright holder.
  18 **      \endlegal
  19 */
  20 /* ========================================================================= */
  21
  22 /* === H E A D E R S ======================================================= */
  23
  24 #ifdef USING_PCH
  25 #       include "pch.h"
  26 #else
  27 #ifdef HAVE_CONFIG_H
  28 #       include <config.h>
  29 #endif
  30
  31 #include "random.h"
  32 #include <cmath>
  33 #include <cstdlib>
  34
  35 #endif
  36
  37 /* === M A C R O S ========================================================= */
  38
  39 /* === G L O B A L S ======================================================= */
  40
  41 /* === P R O C E D U R E S ================================================= */
  42
  43 /* === M E T H O D S ======================================================= */
  44
  45 void
  46 Random::set_seed(int x)
  47 {
  48         seed_=x;
  49         srand(x);
  50         int i;
  51         for(i=0;i<POOL_SIZE;i++)
  52                 pool_[i]=rand();
  53
  54         x_mask=rand()+rand()*RAND_MAX;
  55         y_mask=rand()+rand()*RAND_MAX;
  56         t_mask=rand()+rand()*RAND_MAX;
  57 }
  58
  59 float
  60 Random::operator()(const int salt,const int x,const int y,const int t)const
  61 {
  62         const int salt_hash(pool_[salt&(POOL_SIZE-1)]);
  63
  64         int index(((x^x_mask)+(y^y_mask)*234672+(t^t_mask)*8439573)^salt_hash);
  65
  66         index+=index*(index/POOL_SIZE);
  67
  68         return (float(pool_[index&(POOL_SIZE-1)])/float(RAND_MAX))*2.0f-1.0f;
  69 }
  70
  71 float
  72 Random::operator()(int smooth,int subseed,float xf,float yf,float tf)const
  73 {
  74         int x((int)floor(xf));
  75         int y((int)floor(yf));
  76         int t((int)floor(tf));
  77
  78         switch(smooth)
  79         {
  80         case 4: // cubic
  81                 {
  82                         #define f(j,i,k)        ((*this)(subseed,i,j,k))
  83                         //Using catmull rom interpolation because it doesn't blur at all
  84                         //bezier curve with intermediate ctrl pts: 0.5/3(p(i+1) - p(i-1)) and similar
  85                         float xfa [4], tfa[4];
  86
  87                         //precalculate indices (all clamped) and offset
  88                         const int xa[] = {x-1,x,x+1,x+2};
  89
  90                         const int ya[] = {y-1,y,y+1,y+2};
  91
  92                         const int ta[] = {t-1,t,t+1,t+2};
  93
  94                         const float dx(xf-x);
  95                         const float dy(yf-y);
  96                         const float dt(tf-t);
  97
  98                         //figure polynomials for each point
  99                         const float txf[] =
 100                         {
 101                                 0.5*dx*(dx*(dx*(-1) + 2) - 1),  //-t + 2t^2 -t^3
 102                                 0.5*(dx*(dx*(3*dx - 5)) + 2),   //2 - 5t^2 + 3t^3
 103                                 0.5*dx*(dx*(-3*dx + 4) + 1),    //t + 4t^2 - 3t^3
 104                                 0.5*dx*dx*(dx-1)                                //-t^2 + t^3
 105                         };
 106
 107                         const float tyf[] =
 108                         {
 109                                 0.5*dy*(dy*(dy*(-1) + 2) - 1),  //-t + 2t^2 -t^3
 110                                 0.5*(dy*(dy*(3*dy - 5)) + 2),   //2 - 5t^2 + 3t^3
 111                                 0.5*dy*(dy*(-3*dy + 4) + 1),    //t + 4t^2 - 3t^3
 112                                 0.5*dy*dy*(dy-1)                                //-t^2 + t^3
 113                         };
 114
 115                         const float ttf[] =
 116                         {
 117                                 0.5*dt*(dt*(dt*(-1) + 2) - 1),  //-t + 2t^2 -t^3
 118                                 0.5*(dt*(dt*(3*dt - 5)) + 2),   //2 - 5t^2 + 3t^3
 119                                 0.5*dt*(dt*(-3*dt + 4) + 1),    //t + 4t^2 - 3t^3
 120                                 0.5*dt*dt*(dt-1)                                //-t^2 + t^3
 121                         };
 122
 123                         //evaluate polynomial for each row
 124                         for(int i = 0; i < 4; ++i)
 125                         {
 126                                 for(int j = 0; j < 4; ++j)
 127                                 {
 128                                         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];
 129                                 }
 130                                 xfa[i] = tfa[0]*txf[0] + tfa[1]*txf[1] + tfa[2]*txf[2] + tfa[3]*txf[3];
 131                         }
 132
 133                         //return the cumulative column evaluation
 134                         return xfa[0]*tyf[0] + xfa[1]*tyf[1] + xfa[2]*tyf[2] + xfa[3]*tyf[3];
 135 #undef f
 136                 }
 137                 break;
 138
 139
 140         case 5: // Fast Spline (non-animated)
 141                 {
 142 #define P(x)    (((x)>0)?((x)*(x)*(x)):0.0f)
 143 #define R(x)    ( P(x+2) - 4.0f*P(x+1) + 6.0f*P(x) - 4.0f*P(x-1) )*(1.0f/6.0f)
 144 #define F(i,j)  ((*this)(subseed,i+x,j+y)*(R((i)-a)*R(b-(j))))
 145 #define FT(i,j,k)       ((*this)(subseed,i+x,j+y,k+t)*(R((i)-a)*R(b-(j))*R((k)-c)))
 146 #define Z(i,j) ret+=F(i,j)
 147 #define ZT(i,j,k) ret+=FT(i,j,k)
 148 #define X(i,j)  // placeholder... To make box more symetric
 149 #define XT(i,j,k)       // placeholder... To make box more symetric
 150
 151                 float a(xf-x), b(yf-y);
 152
 153                 // Interpolate
 154                 float ret(F(0,0));
 155                 Z(-1,-1); Z(-1, 0); Z(-1, 1); Z(-1, 2);
 156                 Z( 0,-1); X( 0, 0); Z( 0, 1); Z( 0, 2);
 157                 Z( 1,-1); Z( 1, 0); Z( 1, 1); Z( 1, 2);
 158                 Z( 2,-1); Z( 2, 0); Z( 2, 1); Z( 2, 2);
 159
 160                 return ret;
 161         }
 162
 163         case 3: // Spline (animated)
 164                 {
 165                         float a(xf-x), b(yf-y), c(tf-t);
 166
 167                         // Interpolate
 168                         float ret(FT(0,0,0));
 169                         ZT(-1,-1,-1); ZT(-1, 0,-1); ZT(-1, 1,-1); ZT(-1, 2,-1);
 170                         ZT( 0,-1,-1); ZT( 0, 0,-1); ZT( 0, 1,-1); ZT( 0, 2,-1);
 171                         ZT( 1,-1,-1); ZT( 1, 0,-1); ZT( 1, 1,-1); ZT( 1, 2,-1);
 172                         ZT( 2,-1,-1); ZT( 2, 0,-1); ZT( 2, 1,-1); ZT( 2, 2,-1);
 173
 174                         ZT(-1,-1, 0); ZT(-1, 0, 0); ZT(-1, 1, 0); ZT(-1, 2, 0);
 175                         ZT( 0,-1, 0); XT( 0, 0, 0); ZT( 0, 1, 0); ZT( 0, 2, 0);
 176                         ZT( 1,-1, 0); ZT( 1, 0, 0); ZT( 1, 1, 0); ZT( 1, 2, 0);
 177                         ZT( 2,-1, 0); ZT( 2, 0, 0); ZT( 2, 1, 0); ZT( 2, 2, 0);
 178
 179                         ZT(-1,-1, 1); ZT(-1, 0, 1); ZT(-1, 1, 1); ZT(-1, 2, 1);
 180                         ZT( 0,-1, 1); ZT( 0, 0, 1); ZT( 0, 1, 1); ZT( 0, 2, 1);
 181                         ZT( 1,-1, 1); ZT( 1, 0, 1); ZT( 1, 1, 1); ZT( 1, 2, 1);
 182                         ZT( 2,-1, 1); ZT( 2, 0, 1); ZT( 2, 1, 1); ZT( 2, 2, 1);
 183
 184                         ZT(-1,-1, 2); ZT(-1, 0, 2); ZT(-1, 1, 2); ZT(-1, 2, 2);
 185                         ZT( 0,-1, 2); ZT( 0, 0, 2); ZT( 0, 1, 2); ZT( 0, 2, 2);
 186                         ZT( 1,-1, 2); ZT( 1, 0, 2); ZT( 1, 1, 2); ZT( 1, 2, 2);
 187                         ZT( 2,-1, 2); ZT( 2, 0, 2); ZT( 2, 1, 2); ZT( 2, 2, 2);
 188
 189                         return ret;
 190
 191 /*
 192
 193                         float dx=xf-x;
 194                         float dy=yf-y;
 195                         float dt=tf-t;
 196
 197                         float ret=0;
 198                         int i,j,h;
 199                         for(h=-1;h<=2;h++)
 200                                 for(i=-1;i<=2;i++)
 201                                         for(j=-1;j<=2;j++)
 202                                                 ret+=(*this)(subseed,i+x,j+y,h+t)*(R(i-dx)*R(j-dy)*R(h-dt));
 203                         return ret;
 204 */
 205                 }
 206                 break;
 207 #undef X
 208 #undef Z
 209 #undef F
 210 #undef P
 211 #undef R
 212
 213         case 2: // Cosine
 214         if((float)t==tf)
 215         {
 216                 int x((int)floor(xf));
 217                 int y((int)floor(yf));
 218                 float a=xf-x;
 219                 float b=yf-y;
 220                 a=(1.0f-cos(a*3.1415927))*0.5f;
 221                 b=(1.0f-cos(b*3.1415927))*0.5f;
 222                 float c=1.0-a;
 223                 float d=1.0-b;
 224                 int x2=x+1,y2=y+1;
 225                 return
 226                         (*this)(subseed,x,y,t)*(c*d)+
 227                         (*this)(subseed,x2,y,t)*(a*d)+
 228                         (*this)(subseed,x,y2,t)*(c*b)+
 229                         (*this)(subseed,x2,y2,t)*(a*b);
 230         }
 231         else
 232         {
 233                 float a=xf-x;
 234                 float b=yf-y;
 235                 float c=tf-t;
 236
 237                 a=(1.0f-cos(a*3.1415927))*0.5f;
 238                 b=(1.0f-cos(b*3.1415927))*0.5f;
 239
 240                 // We don't perform this on the time axis, otherwise we won't
 241                 // get smooth motion
 242                 //c=(1.0f-cos(c*3.1415927))*0.5f;
 243
 244                 float d=1.0-a;
 245                 float e=1.0-b;
 246                 float f=1.0-c;
 247
 248                 int x2=x+1,y2=y+1,t2=t+1;
 249
 250                 return
 251                         (*this)(subseed,x,y,t)*(d*e*f)+
 252                         (*this)(subseed,x2,y,t)*(a*e*f)+
 253                         (*this)(subseed,x,y2,t)*(d*b*f)+
 254                         (*this)(subseed,x2,y2,t)*(a*b*f)+
 255                         (*this)(subseed,x,y,t2)*(d*e*c)+
 256                         (*this)(subseed,x2,y,t2)*(a*e*c)+
 257                         (*this)(subseed,x,y2,t2)*(d*b*c)+
 258                         (*this)(subseed,x2,y2,t2)*(a*b*c);
 259         }
 260         case 1: // Linear
 261         if((float)t==tf)
 262         {
 263                 int x((int)floor(xf));
 264                 int y((int)floor(yf));
 265                 float a=xf-x;
 266                 float b=yf-y;
 267                 float c=1.0-a;
 268                 float d=1.0-b;
 269                 int x2=x+1,y2=y+1;
 270                 return
 271                         (*this)(subseed,x,y,t)*(c*d)+
 272                         (*this)(subseed,x2,y,t)*(a*d)+
 273                         (*this)(subseed,x,y2,t)*(c*b)+
 274                         (*this)(subseed,x2,y2,t)*(a*b);
 275         }
 276         else
 277         {
 278
 279                 float a=xf-x;
 280                 float b=yf-y;
 281                 float c=tf-t;
 282
 283                 float d=1.0-a;
 284                 float e=1.0-b;
 285                 float f=1.0-c;
 286
 287                 int x2=x+1,y2=y+1,t2=t+1;
 288
 289                 return
 290                         (*this)(subseed,x,y,t)*(d*e*f)+
 291                         (*this)(subseed,x2,y,t)*(a*e*f)+
 292                         (*this)(subseed,x,y2,t)*(d*b*f)+
 293                         (*this)(subseed,x2,y2,t)*(a*b*f)+
 294                         (*this)(subseed,x,y,t2)*(d*e*c)+
 295                         (*this)(subseed,x2,y,t2)*(a*e*c)+
 296                         (*this)(subseed,x,y2,t2)*(d*b*c)+
 297                         (*this)(subseed,x2,y2,t2)*(a*b*c);
 298         }
 299         default:
 300         case 0:
 301                 return (*this)(subseed,x,y,t);
 302         }
 303 }