ViSP  2.8.0
vpNoise.cpp
1 /****************************************************************************
2  *
3  * \$Id: vpNoise.cpp 4056 2013-01-05 13:04:42Z fspindle \$
4  *
5  * This file is part of the ViSP software.
7  *
8  * This software is free software; you can redistribute it and/or
9  * modify it under the terms of the GNU General Public License
10  * ("GPL") version 2 as published by the Free Software Foundation.
11  * See the file LICENSE.txt at the root directory of this source
13  *
14  * For using ViSP with software that can not be combined with the GNU
17  *
19  *
20  * This software was developed at:
21  * INRIA Rennes - Bretagne Atlantique
22  * Campus Universitaire de Beaulieu
23  * 35042 Rennes Cedex
24  * France
26  *
28  * INRIA at visp@inria.fr
29  *
30  * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
31  * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
32  *
33  *
34  * Description:
35  * Generation of random number with uniform and normal probability density.
36  *
37  * Authors:
38  * Eric Marchand
39  *
40  *****************************************************************************/
41
42 #include <visp/vpNoise.h>
43 #include <math.h>
44
65 inline void
66 vpUniRand::draw0()
67 //minimal standard
68 //Park and Miller random number generator
69 {
70  /*unsigned*/ long k= x/q;//temp value for computing without overflow
71  x = a*(x-k*q)-k*r;
72  if (x < 0) x += m; //compute x without overflow
73 }
74
82 double
84 {
85  const long ntab = 33; //we work on a 32 elements array.
86  //the 33rd one is actually the first value of y.
87  const long modulo = ntab-2;
88
89  static long y = 0;
90  static long T[ntab];
91
92  long j; //index of T
93
94  //step 0
95  if (!y) { //first time
96  for(j = 0; j < ntab; j++) {
97  draw0();
98  T[j]=x;
99  } //compute table T
100  y=T[ntab-1];
101  }
102
103  //step 1
104  j = y & modulo; //compute modulo ntab+1 (the first element is the 0th)
105
106  //step 3
107  y=T[j];
108  double ans = (double)y/normalizer;
109
110  //step 4
111  //generate x(k+i) and set y=x(k+i)
112  draw0();
113
114  //refresh T[j];
115  T[j]=x;
116
117  return ans;
118 }
119
127 double
128 vpGaussRand::gaussianDraw()
129 {
130  double v1, v2, rsq;
131  static bool AlreadyDone = false;
132  static double x2;
133
136  return x2;
137  }
138
139  else {
140
141  do {
142  v1=2*draw1()-1;
143  v2=2*draw1()-1;
144  rsq=v1*v1+v2*v2;
145  } while (rsq >= 1);
146
147  double fac=sqrt(-2*log(rsq)/rsq);
148  x2=v2*fac;