ExampleRandom2D.cpp

/* 
This file is part of the BIAS library (Basic ImageAlgorithmS).

Copyright (C) 2003-2009    (see file CONTACT for details)
  Multimediale Systeme der Informationsverarbeitung
  Institut fuer Informatik
  Christian-Albrechts-Universitaet Kiel


BIAS is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version.

BIAS is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License
along with BIAS; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

/**
   @example ExampleRandom2D.cpp
   @relates Random
   @brief Example2D for random number generator
   @ingroup g_examples
   @author MIP 
*/

#ifdef WIN32
#  include "Base/Common/W32Compat.hh"
#  include <math.h>
#endif

#include <cmath>
#include <cstring>
#include <Base/Math/Random.hh>

using namespace BIAS;
using namespace std;

#define COUNT 9000
/**
   @example ExampleRandom2D.cpp
   @brief Example2D for random number generator
   @ingroup g_examples
   @author MIP 
*/
int main()
{
  Random ran;
  
  ran.Reset();

  // generate random covarinace matrix
  Matrix2x2<double> cov;
  double cov_max=100.0;
  Vector2<double> mean;
  cov.SetZero();

  double sigma1, sigma2;
  sigma1=ran.GetUniformDistributed(0.0, cov_max);
  sigma2=ran.GetUniformDistributed(0.0, cov_max);
  double ang=ran.GetUniformDistributed(0.0, M_PI*2.0);
  //ang=0;
  // covariance is transformed by R * cov * R.Transpose()
  cov[0][0] = cos(ang) * cos(ang) * sigma1 + sin(ang) * sin(ang) * sigma2;
  cov[1][1] = cos(ang) * cos(ang) * sigma2 + sin(ang) * sin(ang) * sigma1;
  cov[0][1] = (sigma1-sigma2)*cos(ang)*sin(ang);
  cov[1][0] = (sigma1-sigma2)*cos(ang)*sin(ang);
  //double a, b, c, d;
  //   a = cov[0][0];
  //   c = cov[1][1];
  //   b = cov[0][1];
  //   d = cov[1][0];
  //  cout << "a "<<a <<"\tb "<<b<<"\tc "<<c<<"\td "<<d<<endl;
  //   cout << sigma1*sigma1+sigma1*(-c-a)+(a*c-b*d) << endl;
  //   cout << sigma2*sigma2+sigma2*(-c-a)+(a*c-b*d) << endl;
  //   cout << (a-sigma1)*(c-sigma1)-b*d << endl;
  //   cout << (a-sigma2)*(c-sigma2)-b*d << endl;

  mean[0] = ran.GetUniformDistributed(-cov_max, cov_max);
  mean[1] = ran.GetUniformDistributed(-cov_max, cov_max);
  //mean.SetZero();
  cout << "input mean: "<<mean<<endl;
  //cout << "sigma1: "<<sigma1<<"\tsigma2: "<<sigma2<<"\tangle: "<<ang<<endl;
  cout << "input covariance: "<<cov<<endl;

  vector<Vector2<double> > points;

  if (ran.GetNormalDistributed(mean, cov, COUNT, points)!=0){
    BIASERR("error");
    return -1;
  }

  // now calculate
  double dx, dy;
  Matrix2x2<double> ocov;
  Vector2<double> omean;

  // zero mean and var
  omean.SetZero();
  ocov.SetZero();
  

  // calculate mean
  for (register unsigned int i=0; i<COUNT; i++){
    //cout << points[i] << endl;
    omean[0]+=points[i][0];
    omean[1]+=points[i][1];
  }
  
  omean[0]/=COUNT;
  omean[1]/=COUNT;
  cout << "output mean:" << omean << endl;

  // calculate covariance matrix
  for (register unsigned int i=0; i<COUNT; i++){
    dx = points[i][0]-omean[0];
    dy = points[i][1]-omean[1];
    ocov[0][0] += dx*dx;
    ocov[1][1] += dy*dy;
    ocov[0][1] += dx*dy;
  }
  ocov[1][0]=ocov[0][1];
  ocov/=(double)(COUNT-1);

  // output of the covariance matrix
  cout << "output covariance:" << ocov << endl;

  return 0;
}