Bingham_Watson_approx.cc

/*  Directional Statistics Functions

    Bingham and Watson Distributions and functions to approximate their normalizing constant
    
    Stam Sotiropoulos  - FMRIB Image Analysis Group
 
    Copyright (C) 2011 University of Oxford  */

/*  CCOPYRIGHT  */

#include "Bingham_Watson_approx.h"


  //Cubic root
  float croot(const float& x){
    float res;
    if (x>=0) 
      res=pow(x,INV3);
    else
      res=-(pow(-x,INV3));
    return res;
  }


  //Saddlepoint approximation of confluent hypergeometric function of a matrix argument B (Kume & Wood, 2005)
  //Vector x has the eigenvalues of B.
  float hyp_Sapprox(ColumnVector &x){
    float c1;

    //SortDescending(x);         //Not needed??
    if (x(1)==0 && x(2)==0 && x(3)==0)
      c1=1;
    else{
      float t=find_t(x);
      float T,R=1, K2=0, K3=0, K4=0;
      for (int i=1; i<=3; i++){
        R=R/sqrt(-x(i)-t);
	K2=K2+1.0/(2*pow(x(i)+t,2.0));
	K3=K3-1.0/pow(x(i)+t,3.0);
	K4=K4+3.0/pow(x(i)+t,4.0);
      }
      T=K4/(8*K2*K2)-5*K3*K3/(24*K2*K2*K2);  
      //c1=sqrt(2.0/K2)*pi*R*exp(-t);
      //float c3=c1*exp(T);
      //c1=c3/(4*pi);
      c1=0.25*sqrt(2.0/K2)*R*exp(-t+T);
    }
    return c1;
  }


  //Saddlepoint approximation of confluent hypergeometric function of a matrix argument, with its eigenvalues being l1,l1,l2 or l1,l2,l2 with l1!=l2.
  //Vector x has the three eigenvalues. This function can be also used to approximate a confluent hypergeometric function of a scalar argument k
  //by providing x=[k 0 0].
  float hyp_Sapprox_twoequal(ColumnVector& x){
    float c1, R, t, Bm, Bp, K2, K3, K4,T,q;

    //SortDescending(x);  //Not needed??
    if (x(1)==x(2)){
      q=1;
      x(2)=x(3);
    }
    else 
      q=2;

    R=sqrt(4*(x(1)-x(2))*(x(1)-x(2))+9+4*(2*q-3)*(x(2)-x(1)));
    t=(-2*(x(1)+x(2))-3-R)/4.0;
    Bm=1/(R+3-2*(x(2)-x(1)));
    Bp=1/(R+3+2*(x(2)-x(1)));
    K2=1.0/(2*pow(x(1)+t,2.0))+1.0/pow(x(2)+t,2.0);