Index: /GTP/trunk/Lib/Vis/Preprocessing/src/Vector2.cpp
===================================================================
--- /GTP/trunk/Lib/Vis/Preprocessing/src/Vector2.cpp	(revision 1975)
+++ /GTP/trunk/Lib/Vis/Preprocessing/src/Vector2.cpp	(revision 1975)
@@ -0,0 +1,145 @@
+// ===================================================================
+// $Id: vector2.cpp,v 1.3 2005/11/30 15:46:50 havran Exp $
+//
+// vector2.cpp
+//     CVector2D class implements an 2D vector
+//
+// Licence: the use and distribution of this file is severely limited, please
+// see the file 'doc/Licence.txt'. Any non-authorized use can be prosecuted under
+// International Law. For further questions, please, e-mail to VHavran@seznam.cz
+// or mail to Vlastimil Havran, Pohodli 27, 57001 Litomysl, the Czech Republic.
+// REPLACEMENT_STRING
+//
+// Initial coding by Jiri Bittner
+
+// GOLEM headers
+#include "Vector2.h"
+
+namespace GtpVisibilityPreprocessor {
+
+// Overload << operator for C++-style output
+ostream&
+operator<< (ostream &s, const Vector2 &A)
+{
+  return s << "(" << A.x() << ", " << A.y() << ")";
+}
+
+// Overload >> operator for C++-style input
+istream&
+operator>> (istream &s, Vector2 &A)
+{
+  char a;
+  // read "(x, y)"
+  return s >> a >> A.x() >> a >> A.y() >> a;
+}
+
+
+// Precompute which.  which is 0 if the normal is dominant
+// in the X direction, 1 if the Y direction.
+int
+Vector2::DominantAxis()
+{
+  return (fabs(x()) > fabs(y())) ?  0 : 1;
+}
+
+void
+Vector2::ExtractVerts(float *p, int dominant) const
+{
+  if (dominant == 0)
+    *p = y();
+  else
+    *p = x();
+}
+
+float
+Vector2::Normalize()
+{
+  float s = Size();
+
+  if (s != 0.0) {
+    xx /= s;
+    yy /= s;
+  }
+  return s;
+}
+
+Vector2
+Vector2::operator-(const Vector2 &v) const
+{
+  Vector2 u;
+  u.SetX(x()-v.x());
+  u.SetY(y()-v.y());
+  return u;
+}
+
+Vector2
+Vector2::operator+(const Vector2 &v) const
+{
+  Vector2 u;
+  u.SetX(x()+v.x());
+  u.SetY(y()+v.y());
+  return u;
+}
+
+Vector2
+Vector2::operator*(const float t) const
+{
+  Vector2 u;
+  u.SetX(t*x());
+  u.SetY(t*y());
+  return u;
+}
+
+int
+Vector2::Equal(const Vector2 &u, float trash) const
+{
+ return ( (fabs(x() - u.x()) < trash) && (fabs(y() - u.y()) < trash) ); 
+}
+
+
+// get the angle between two vectors in range 0 - PI
+
+float
+Vector2::Angle(const Vector2 &v) const
+{
+  float cosine;
+  cosine = DotProd(*this, v) / (Size() * v.Size());
+  return acos(cosine);
+}
+
+float
+Vector2::Cosine(const Vector2 &v) const
+{
+  return DotProd(*this, v)  /  (Size() * v.Size());
+}
+
+// cosine assuming that this vector is normalized
+float
+Vector2::CosineN(const Vector2 &v) const
+{
+  return DotProd(*this, v) / v.Size();
+}
+
+// if a given vector is smaller in one of coordinates than this, update
+Vector2&
+Vector2::UpdateMin(const Vector2 &v)
+{
+  if (x() > v.x())
+    SetX(v.x());
+  if (y() > v.y())
+    SetY(v.y());
+  return *this;
+}
+
+Vector2&
+Vector2::UpdateMax(const Vector2 &v)
+{
+  if ( x() < v.x() )
+    SetX(v.x());
+  if ( y() < v.y() )
+    SetY(v.y());
+  return *this;
+}
+
+
+}
Index: /GTP/trunk/Lib/Vis/Preprocessing/src/RndGauss.h
===================================================================
--- /GTP/trunk/Lib/Vis/Preprocessing/src/RndGauss.h	(revision 1975)
+++ /GTP/trunk/Lib/Vis/Preprocessing/src/RndGauss.h	(revision 1975)
@@ -0,0 +1,141 @@
+// ===================================================================
+// $Id:$
+//
+// rndgauss.h:
+//     A set of functions to generate random numbers according to
+//     bivariate normal (=Gaussian) distribution + special distribution
+//     on disk combining it with bivariate normal distribution
+//
+// Initial coding by Vlastimil Havran, January 2007
+
+#ifndef __RNDGAUSS_H__
+#define __RNDGAUSS_H__
+
+//#include "configh.h"
+//#include "basrnd.h"
+//#include "errorg.h"
+#include <assert.h>
+#include "Vector2.h"
+//#include "sphcoord.h"
+
+namespace GtpVisibilityPreprocessor {
+
+// Function: GaussianOn2D
+// Here we generate the Gaussian distribution on the 2D domain
+// (=bivariate normal distribution) given the uniformly distributed random
+// numbers in the range [0-1)x[0-1) (bivariate uniform distribution).
+// We specify also standard deviation (sigma) of the normal distribution.
+// The result has a mean value <0.0f, 0.0f>. We use Box-Muller transform
+// in polar form.
+// Please note that in the distance from mean value this is true:
+// the distance = 1 x sigma ... we have 68.25% of generated points
+//              = 2 x sigma ... we have 95.45% of generated points
+//              = 3 x sigma ... we have 99.73% of generated points
+//              = 4 x sigma ... we have 99.99% of generated points
+
+inline void
+GaussianOn2D(const Vector2 &input, // input (RND in interval [0-1)x[0-1))
+			 float sigma,            // standard deviation of gaussian distribution
+			 Vector2 &outputVec)   // resulting vector according to gaussian distr.
+{
+  // Handling zero
+  float xf = input.xx;
+  if (xf < 1e-20f)
+    xf = 1e-20f;
+  assert( xf <= 1.0f);
+  assert( (input.yy >= 0.f) && (input.yy <= 1.0f));
+
+  // Here we have to use natural logarithm function !
+  float t = -2.0f * log(xf);
+  assert(t > 0.f);
+  t = sigma * sqrt(t);
+  float phi = 2.0f * M_PI * input.yy;
+  outputVec.xx = t * cos(phi);
+  outputVec.yy = t * sin(phi);
+
+  return;
+}
+
+// Function: PolarGaussianOnDisk
+// Here we generate the Gaussian distribution on the 2D domain on the part
+// of 2D disk. The input of the function is tuple of random
+// numbers in the range [0-1)x[0-1) (bivariate uniform distribution).
+// The result of this function 'outputVec' (in general in the range
+// [-infinity,+infinity]x[-infinity,+infinity], so the point generated
+// according to proposed distribution. The algorithm does not use rejection
+// sampling. We have to specify the radius of 2D circle
+// where the distribution mean value is zero for one variable !!!!
+// Sigma is specified independently, but if you want to have small probability
+// of generating the original point, then consider:
+// the radius = 1 x sigma ... 68.25% points in the distance from sigma from the circle
+//            = 2 x sigma ... 95.45% points in the distance from sigma from the circle
+//            = 3 x sigma ... 99.73% points in the distance from sigma from the circle
+//            = 4 x sigma ... 99.99% points in the distance from sigma from the circle
+// Recommendation: use sigma = radius/3.3f or so (according to visualization)
+
+// This is only helper function
+inline float
+ComputeSigmaFromRadius(float radius)
+{
+  // For standard deviation = 3 we have achieve probability 0.4% that the proximity
+  // of the center of the circle will be covered again by a generated point.
+  return radius / 3.3f;
+}
+
+inline void
+PolarGaussianOnDisk(const Vector3 &input, // input (RND in interval [0-1)x[0-1))
+					float sigma,    // standard deviation of gaussian distribution
+					float radius,   // standard deviation of gaussian distribution
+					Vector2 &outputVec) // result
+{
+  // Handling zero
+  float xf = input.x;
+  if (xf < 1e-20f)
+    xf = 1e-20f;
+  assert( xf <= 1.0f);
+  assert( (input.y >= 0.f) && (input.y <= 1.0f));
+
+  // Here we have to use natural logarithm function !
+  float t = -2.0f * log(xf);
+  assert(t > 0.f);
+  t = sigma * sqrt(t);
+  float phi = 2.0f * M_PI * input.y;
+
+#if 1
+  // This is what Jiri Bittner wanted originally,
+  // another idea how to implement this distribution
+  // First, we generate point as for standard polar
+  // inmplementation
+  outputVec.xx = t * cos(phi);
+  outputVec.yy = t * sin(phi);
+
+  // Then we generate a point on 2D circle at the distance
+  // of radius
+  float r = input.z; //RND::rand01f();
+  float phi2 = 2.0f * M_PI * r;
+  Vector2 c(cos(phi2),sin(phi2));
+  outputVec.xx += radius * c.xx;
+  outputVec.yy += radius * c.yy;
+#endif
+#if 0
+  // This is what Jiri Bittner wanted originally
+  // however, I cannot make it running correctly !!!
+  // Do not use it at the moment !!! VH 11/1/2007
+  float r = RND::rand01f();  
+  if (r < 0.5f) {
+    outputVec.xx = (radius + t) * cos(phi);
+    outputVec.yy = (radius + t) * sin(phi);
+  }
+  else {
+    // going inside the cirle
+    outputVec.xx = (radius - t) * cos(phi);
+    outputVec.yy = (radius - t) * sin(phi);    
+  }
+#endif
+  
+  return;
+}
+
+}
+
+#endif // __RND_H__
Index: /GTP/trunk/Lib/Vis/Preprocessing/src/Vector2.h
===================================================================
--- /GTP/trunk/Lib/Vis/Preprocessing/src/Vector2.h	(revision 1975)
+++ /GTP/trunk/Lib/Vis/Preprocessing/src/Vector2.h	(revision 1975)
@@ -0,0 +1,202 @@
+ // ===================================================================
+// $Id: vector2.h,v 1.5 2005/11/02 12:58:48 xkrivanj Exp $
+//
+// vector2.h
+//     Header file for CVector2D class - implements 2-dimensional vector
+//
+// Class: CVector2D
+//
+// Licence: the use and distribution of this file is severely limited, please
+// see the file 'doc/Licence.txt'. Any non-authorized use can be prosecuted under
+// International Law. For further questions, please, e-mail to VHavran@seznam.cz
+// or mail to Vlastimil Havran, Pohodli 27, 57001 Litomysl, the Czech Republic.
+// REPLACEMENT_STRING
+//
+// Initial coding by Jiri Bittner
+
+#ifndef __VECTOR2_H__
+#define __VECTOR2_H__
+
+// #include "basmacr.h"
+//#include "basmath.h"
+#include "common.h"
+#include <iostream>
+using namespace std;
+
+namespace GtpVisibilityPreprocessor {
+
+// ---------------------------------------------------------------------------
+// This class describes two-dimensional vector.
+// ---------------------------------------------------------------------------
+class Vector2
+{
+public:
+  float xx, yy; // coordinates of the vector
+
+  // constructors
+  Vector2(float xs, float ys):xx(xs), yy(ys) {}
+  // default constructor
+  Vector2() {} 
+  
+  friend ostream& operator<< (ostream &s, const Vector2 &A);
+  friend istream& operator>> (istream &s, Vector2 &A);
+
+  const   float&   x() const { return xx; }
+  const   float&   y() const { return yy; }
+  float&   x() {return xx;}
+  float&   y() {return yy;}
+  
+  // Functions to get at the vector components
+  float& operator[] (int inx) {
+    if (inx == 0) return xx;
+    else return yy;
+  }
+
+  const float& operator[] (int inx) const {
+    if (inx == 0) return xx;
+    else return yy;
+  }
+
+  // returns x-coordinate for which==0, x-coordinate for which==1
+  void     ExtractVerts(float *p, int which) const;
+
+  float&   SetX(const float x) {return xx = x;}
+  float&   SetY(const float y) {return yy = y;}
+
+  Vector2& operator=(Vector2 const &v) {
+    xx = v.x(); yy = v.y(); return *this;
+  }
+
+  void   SetValue(float xs, float ys ) { xx = xs; yy = ys; }
+
+  float  Size() const { return (float) sqrt(xx*xx + yy*yy);}
+
+  // returns the squared magnitude of a vector
+  friend inline float SqrMagnitude(const Vector2 &v);
+
+  // returns the squared distance between two vectors
+  friend inline float SqrDistance(const Vector2 &v1, const Vector2 &v2);
+  
+  // normalize the vector to the size=1.0
+  float    Normalize();
+  
+  // Assignment operators
+  Vector2& operator+= (const Vector2 &a);
+  Vector2& operator-= (const Vector2 &a);
+  Vector2& operator*= (const Vector2 &a);
+  Vector2& operator*= (float a);
+  Vector2& operator/= (float a);
+
+  Vector2 operator-(const Vector2 &) const;
+  Vector2 operator+(const Vector2 &) const;
+  Vector2 operator*(const float t) const;
+
+  //  Vector2 operator*(const Matrix3C &mat) const;
+  Vector2 operator-() const { return Vector2(-x(),-y()); }
+
+  int  operator==(const Vector2 &u) const {return Equal(u,Limits::Small);}
+  int  operator!=(const Vector2 &u) const {return !operator == (u);}
+
+  int  Equal(const Vector2 &u,float trash) const;
+
+  // dot product of the two vectors
+  friend inline double DotProd(const Vector2 &u,
+			       const Vector2 &v) {
+    return ( u.x() * v.x()  +  u.y() * v.y() );
+  }
+
+  // the angle between two vectors $\in <0, PI>$
+  float  Angle(const Vector2 &v) const;
+  // cosine of the angle between the two vectors $\in <-1,1>$
+  float  Cosine(const Vector2 &v) const;
+  // supposes this vector is normalized
+  float  CosineN(const Vector2 &v) const;
+
+  // checks if the this vector is not opposite to a given vector
+  int    IsOpposite(const Vector2 &v) const {
+    return (Abs(DotProd(*this,v) + Size() * v.Size()) < Limits::Small);
+  }
+
+  // computes the distance between this and a given vector
+  float  Distance(const Vector2 &v) const { 
+    return  (float) sqrt( sqr(xx-v.x()) + sqr(yy - v.y()));
+  }
+
+  // computes the squared distance between this and a given vector 
+  float  SqrDistance(const Vector2 &v) const { 
+    return  sqr(x()-v.x())+sqr(y()-v.y());
+  }
+
+  // if a given vector has a smaller(larger) coordinate, then this is updated
+  Vector2&  UpdateMin(const Vector2 &v);
+  Vector2&  UpdateMax(const Vector2 &v);
+
+  // checks if both coordinates are smaller than a given one
+  int operator<=(Vector2 &v) {
+    return x() < (v.x() + Limits::Small)
+      && y() < (v.y()+Limits::Small);
+  }
+
+  // checks if both coordinates are larger than a given one
+  int operator>=(Vector2 &v) {
+    return x() > (v.x()-Limits::Small) &&
+      y() > (v.y()-Limits::Small);
+  }
+
+  // returns which coordinates has larger size
+  int DominantAxis();
+};
+
+const Vector2 ZeroVector2(0,0);
+
+inline Vector2&
+Vector2::operator+= (const Vector2 &a)
+{
+  xx += a.xx;  yy += a.yy;
+  return *this;
+}
+
+inline Vector2&
+Vector2::operator-= (const Vector2 &a)
+{
+  xx -= a.xx;  yy -= a.yy;
+  return *this;
+}
+
+inline Vector2&
+Vector2::operator*= (float a)
+{
+  xx *= a;  yy *= a;
+  return *this;
+}
+
+inline Vector2&
+Vector2::operator/= (float a)
+{
+  xx /= a;  yy /= a;
+  return *this;
+}
+
+inline Vector2&
+Vector2::operator*= (const Vector2 &a)
+{
+  xx *= a.xx;  yy *= a.yy;
+  return *this;
+}
+
+inline float
+SqrMagnitude(const Vector2 &v)
+{
+  return v.x() * v.x() + v.y() * v.y();
+}
+
+inline float
+SqrDistance(const Vector2 &v1, const Vector2 &v2)
+{
+  return sqr(v1.xx - v2.xx) + sqr(v1.yy - v2.yy);
+}
+
+}
+
+#endif // __VECTOR2_H__
+
