ViSP  2.9.0
vpNoise.cpp
1 /****************************************************************************
2  *
3  * $Id: vpNoise.cpp 4574 2014-01-09 08:48:51Z fspindle $
4  *
5  * This file is part of the ViSP software.
6  * Copyright (C) 2005 - 2014 by INRIA. All rights reserved.
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
12  * distribution for additional information about the GNU GPL.
13  *
14  * For using ViSP with software that can not be combined with the GNU
15  * GPL, please contact INRIA about acquiring a ViSP Professional
16  * Edition License.
17  *
18  * See http://www.irisa.fr/lagadic/visp/visp.html for more information.
19  *
20  * This software was developed at:
21  * INRIA Rennes - Bretagne Atlantique
22  * Campus Universitaire de Beaulieu
23  * 35042 Rennes Cedex
24  * France
25  * http://www.irisa.fr/lagadic
26  *
27  * If you have questions regarding the use of this file, please contact
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 
134  if (AlreadyDone) {
135  AlreadyDone=false;
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;
149  AlreadyDone=true;
150  return v1*fac;
151  }
152 }
153 
double draw1()
Bays-Durham Shuffling of Park-Miller generator.
Definition: vpNoise.cpp:83
long x
Definition: vpNoise.h:83