ViSP
Main Page
Related Pages
Modules
Classes
Examples
All
Classes
Functions
Variables
Enumerations
Enumerator
Friends
Groups
Pages
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.
6
* Copyright (C) 2005 - 2013 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
83
vpUniRand::draw1
()
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
src
math
misc
vpNoise.cpp
Generated on Tue Sep 17 2013 00:21:36 for ViSP by
1.8.4