1 | // ===================================================================
|
---|
2 | // $Id:$
|
---|
3 | //
|
---|
4 | // rndgauss.h:
|
---|
5 | // A set of functions to generate random numbers according to
|
---|
6 | // bivariate normal (=Gaussian) distribution + special distribution
|
---|
7 | // on disk combining it with bivariate normal distribution
|
---|
8 | //
|
---|
9 | // Initial coding by Vlastimil Havran, January 2007
|
---|
10 |
|
---|
11 | #ifndef __RNDGAUSS_H__
|
---|
12 | #define __RNDGAUSS_H__
|
---|
13 |
|
---|
14 | #include "SampleGenerator.h"
|
---|
15 | #include "Vector3.h"
|
---|
16 |
|
---|
17 | #include <assert.h>
|
---|
18 |
|
---|
19 |
|
---|
20 | namespace CHCDemoEngine
|
---|
21 | {
|
---|
22 | // Function: GaussianOn2D
|
---|
23 | // Here we generate the Gaussian distribution on the 2D domain
|
---|
24 | // (=bivariate normal distribution) given the uniformly distributed random
|
---|
25 | // numbers in the range [0-1)x[0-1) (bivariate uniform distribution).
|
---|
26 | // We specify also standard deviation (sigma) of the normal distribution.
|
---|
27 | // The result has a mean value <0.0f, 0.0f>. We use Box-Muller transform
|
---|
28 | // in polar form.
|
---|
29 | // Please note that in the distance from mean value this is true:
|
---|
30 | // the distance = 1 x sigma ... we have 68.25% of generated points
|
---|
31 | // = 2 x sigma ... we have 95.45% of generated points
|
---|
32 | // = 3 x sigma ... we have 99.73% of generated points
|
---|
33 | // = 4 x sigma ... we have 99.99% of generated points
|
---|
34 |
|
---|
35 | inline void
|
---|
36 | GaussianOn2D(const Sample2 &input, // input (RND in interval [0-1)x[0-1))
|
---|
37 | float sigma, // standard deviation of gaussian distribution
|
---|
38 | Sample2 &outputVec) // resulting vector according to gaussian distr.
|
---|
39 | {
|
---|
40 | // Handling zero
|
---|
41 | float xf = input.x;
|
---|
42 | if (xf < 1e-20f)
|
---|
43 | xf = 1e-20f;
|
---|
44 | assert( xf <= 1.0f);
|
---|
45 | assert( (input.y >= 0.f) && (input.y <= 1.0f));
|
---|
46 |
|
---|
47 | // Here we have to use natural logarithm function !
|
---|
48 | float t = -2.0f * log(xf);
|
---|
49 | assert(t > 0.f);
|
---|
50 | t = sigma * sqrt(t);
|
---|
51 | float phi = 2.0f * M_PI * input.y;
|
---|
52 | outputVec.x = t * cos(phi);
|
---|
53 | outputVec.y = t * sin(phi);
|
---|
54 |
|
---|
55 | return;
|
---|
56 | }
|
---|
57 |
|
---|
58 | // Function: PolarGaussianOnDisk
|
---|
59 | // Here we generate the Gaussian distribution on the 2D domain on the part
|
---|
60 | // of 2D disk. The input of the function is tuple of random
|
---|
61 | // numbers in the range [0-1)x[0-1) (bivariate uniform distribution).
|
---|
62 | // The result of this function 'outputVec' (in general in the range
|
---|
63 | // [-infinity,+infinity]x[-infinity,+infinity], so the point generated
|
---|
64 | // according to proposed distribution. The algorithm does not use rejection
|
---|
65 | // sampling. We have to specify the radius of 2D circle
|
---|
66 | // where the distribution mean value is zero for one variable !!!!
|
---|
67 | // Sigma is specified independently, but if you want to have small probability
|
---|
68 | // of generating the original point, then consider:
|
---|
69 | // the radius = 1 x sigma ... 68.25% points in the distance from sigma from the circle
|
---|
70 | // = 2 x sigma ... 95.45% points in the distance from sigma from the circle
|
---|
71 | // = 3 x sigma ... 99.73% points in the distance from sigma from the circle
|
---|
72 | // = 4 x sigma ... 99.99% points in the distance from sigma from the circle
|
---|
73 | // Recommendation: use sigma = radius/3.3f or so (according to visualization)
|
---|
74 |
|
---|
75 | // This is only helper function
|
---|
76 | inline float
|
---|
77 | ComputeSigmaFromRadius(float radius)
|
---|
78 | {
|
---|
79 | // For standard deviation = 3 we have achieve probability 0.4% that the proximity
|
---|
80 | // of the center of the circle will be covered again by a generated point.
|
---|
81 | return radius / 3.3f;
|
---|
82 | }
|
---|
83 |
|
---|
84 | inline void
|
---|
85 | PolarGaussianOnDisk(const Vector3 &input, // input (RND in interval [0-1)x[0-1))
|
---|
86 | float sigma, // standard deviation of gaussian distribution
|
---|
87 | float radius, // standard deviation of gaussian distribution
|
---|
88 | Sample2 &outputVec) // result
|
---|
89 | {
|
---|
90 | // Handling zero
|
---|
91 | float xf = input.x;
|
---|
92 | if (xf < 1e-20f)
|
---|
93 | xf = 1e-20f;
|
---|
94 | assert( xf <= 1.0f);
|
---|
95 | assert( (input.y >= 0.f) && (input.y <= 1.0f));
|
---|
96 |
|
---|
97 | // Here we have to use natural logarithm function !
|
---|
98 | float t = -2.0f * log(xf);
|
---|
99 | assert(t > 0.f);
|
---|
100 | t = sigma * sqrt(t);
|
---|
101 | float phi = 2.0f * M_PI * input.y;
|
---|
102 |
|
---|
103 | #if 1
|
---|
104 | // This is what Jiri Bittner wanted originally,
|
---|
105 | // another idea how to implement this distribution
|
---|
106 | // First, we generate point as for standard polar
|
---|
107 | // inmplementation
|
---|
108 | outputVec.x = t * cos(phi);
|
---|
109 | outputVec.y = t * sin(phi);
|
---|
110 |
|
---|
111 | // Then we generate a point on 2D circle at the distance
|
---|
112 | // of radius
|
---|
113 | float r = input.z; //RND::rand01f();
|
---|
114 | float phi2 = 2.0f * M_PI * r;
|
---|
115 |
|
---|
116 | Sample2 c(cos(phi2),sin(phi2));
|
---|
117 |
|
---|
118 | outputVec.x += radius * c.x;
|
---|
119 | outputVec.y += radius * c.y;
|
---|
120 |
|
---|
121 | #endif
|
---|
122 | #if 0
|
---|
123 | // This is what Jiri Bittner wanted originally
|
---|
124 | // however, I cannot make it running correctly !!!
|
---|
125 | // Do not use it at the moment !!! VH 11/1/2007
|
---|
126 | float r = RND::rand01f();
|
---|
127 | if (r < 0.5f) {
|
---|
128 | outputVec.x = (radius + t) * cos(phi);
|
---|
129 | outputVec.y = (radius + t) * sin(phi);
|
---|
130 | }
|
---|
131 | else {
|
---|
132 | // going inside the cirle
|
---|
133 | outputVec.x = (radius - t) * cos(phi);
|
---|
134 | outputVec.y = (radius - t) * sin(phi);
|
---|
135 | }
|
---|
136 | #endif
|
---|
137 |
|
---|
138 | return;
|
---|
139 | }
|
---|
140 |
|
---|
141 | }
|
---|
142 | #endif // __RND_H__
|
---|