source: GTP/trunk/Lib/Vis/Preprocessing/src/RndGauss.h @ 1975

Revision 1975, 4.8 KB checked in by bittner, 18 years ago (diff)

gauss sampling

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