synfig-core/tags/stable/src/modules/mod_noise/random.cpp

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