freAsymmetricVariateGenerator.cxx

Go to the documentation of this file.

/*=========================================================================

  Program:   F.R.E.E. - flexible registration evaluation engine
  Version:   v.1.0.0
  Date:      $Date: 2006/09/01 12:00:00 $
  Module:    $RCSfile: freAsymmetricVariateGenerator.cxx,v $


  Copyright (c) 2007 Ralf o Floca (Department of Medical Informatics,
  Institute for Medical Biometry and Informatics, University of Heidelberg,
  Germany). All rights reserved.
  See FREECopyright.txt or http://www.mi.med.uni-hd.de/free/copyright.htm
  for details.

     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/
#include "freAsymmetricVariateGenerator.h"
#include "vnl/vnl_sample.h"
#include "vnl/vnl_math.h"

namespace itk {
  namespace Statistics {

    void
      AsymmetricVariateGenerator::
      PrintSelf(std::ostream& os, Indent indent) const
    {
      Superclass::PrintSelf(os, indent) ;

      // Print state vector contents
      os << indent << "Variance: " << m_Variance << std::endl;
      os << indent << "Skewness: " << m_Skewness << std::endl;
      os << indent << "Gamma: " << m_Gamma << std::endl;
    }

      AsymmetricVariateGenerator::
      AsymmetricVariateGenerator()
    {
      Initialize();
    };

    void
      AsymmetricVariateGenerator::
      Initialize( const int randomSeed )
    {
      vnl_sample_reseed(randomSeed);
    };

    void
      AsymmetricVariateGenerator::
      Initialize()
    {
      vnl_sample_reseed();
    };  

    double
      AsymmetricVariateGenerator::
      GetVariate()
    {
      /*implementation of the function asym (see Reference)*/

      double g = vnl_sample_uniform(0,1.0);
      double t = pow(1+m_Skewness , m_Gamma/2);
      double gTreshold = 1/(1+t);

      double result = 0.0;

      if (g<gTreshold)
      {
        result = sqrt(2*pow(m_Variance,m_Gamma)) * Erf((t*g)+g-1);
      }
      else
      {
        double h = pow(1+m_Skewness,-1*(m_Gamma/2));
        result = sqrt(2*pow((1+m_Skewness)*m_Variance , m_Gamma)) * Erf((h*g)-h+g);
      };
      return result;
    };

    double
      AsymmetricVariateGenerator::
      operator()()
    {
      return this->GetVariate();
    }; 

    double
      AsymmetricVariateGenerator::
      Erf(double x)
    {
      double result = (1/sqrt(2.0)) * INorm((x+1)/2);
      return result;
    };

    double
      AsymmetricVariateGenerator::
      Freq(double x)
    {
      double v = abs(x);
      if (v<0.5)
      {
        double h = (3.6767877*v - 0.097970465*pow(v,3))/(3.2584593 + pow(v,2));
        if (x>0) return 0.5*h + 0.5;
        else return 0.5 - 0.5*h;
      }
      else if (v<4.0)
      {
        //precalculate, because it is used twice
        double v_2 = v*v; 
        double v_3 = v_2*v;
        double v_4 = v_3*v;

        double ap = 7.37389 + 6.86502*v + 3.0318*v_2 + 0.56317*v_3 + 0.0000431878*v_4;
        double aq = 7.37396 + 15.1849*v + 12.7955*v_2 + 5.35422*v_3 + v_4;
        double h = 1 - exp(-1*v*v) * (ap/aq);
        if (x>0) return 0.5*h + 0.5;
        else return 0.5 - 0.5*h;
      }
      else
      {
        double h = 1 - exp(-1*v*v) * ((-0.96821 + 0.439821*pow(v,2) + 0.248761*pow(v,4)) / (pow(v,3) + 0.440917*pow(v,5)));
        if (x>0) return 0.5*h + 0.5;
        else return 0.5 - 0.5*h;
      }
    };

    double
      AsymmetricVariateGenerator::
      INorm(double z)
    {
      if ((z<=0) || (z>=1)) itkExceptionMacro(<<"Invalid input value z. Value must be in interval (0, 1). z = "<<z);

      if (z==0.5) return 0;

      double x = z;

      if (z>0.5) x = 1-z;

      x = sqrt(-2*log(x));
      //seems to be a typo in the original reference, there the term "450.636*x²" is used.
      //As far as I know the original algorithm by Hill and Davis uses 450.636*x.
      x = x - ((1271.06 + 450.636*x + 7.45551*pow(x,2)) / (500.756 + 750.365*x + 110.421*pow(x,2) + pow(x,3)));

      if (z<0,5) x = -1 * x;

      return x + 2.50662827*(z-Freq(x))*exp(0.5*x*x);
    }
    ;
  }

}

---