Initial revision

base2z.cc

+/*  base2z.cc
+
+    Mark Woolrich & Mark Jenkinson, FMRIB Image Analysis Group
+
+    Copyright (C) 1999-2000 University of Oxford  */
+
+/*  CCOPYRIGHT  */
+
+#include <cmath>
+#include "base2z.h"
+#include "libprob.h"
+
+namespace MISCMATHS {
+
+  Base2z* Base2z::base2z = NULL;
+
+  float Base2z::logbeta(float v, float w)
+  {
+    return MISCMATHS::lgam(v)+MISCMATHS::lgam(w)-MISCMATHS::lgam(v+w);
+  }
+
+  float Base2z::logp2largez(float logp) 
+  {
+    // Large Z extrapolation routine for converting log(p) to Z values
+    //  written by Mark Jenkinson, March 2000
+    //
+    // Equations were derived by using integration by parts and give the
+    //  following formulae:
+    //   Z to log(p)
+    //       log(p) = -1/2*z*z - 1/2*log(2*pi) - log(z) 
+    //                + log(1 - 1/(z*z) + 3/(z*z*z*z))
+    // this equation is then solved by the recursion:
+    //   z_0 = sqrt(2*(-log(p) - 1/2*log(2*pi)))
+    //   z_{n+1} = sqrt(2*(-log(p) - 1/2*log(2*pi) - log(z_n)  
+    //             + log(1 - 1/(zn*zn) + 3/(zn*zn*zn*zn)) ))
+    // In practice this recursion is quite accurate in 3 to 5 iterations
+    // The equation is accurate to 1 part in 10^3 for Z>3.12  (3 iterations)
+
+
+    static const float pi = 3.141592653590;
+    static const float log2pi = log(2*pi);
+    
+    float z0, zn;
+    // iteratively solve for z given log p
+    float b = -2*logp - log2pi; 
+    z0 = sqrt(b);
+    zn = z0;
+    for (int m=1; m<=3; m++) {
+      // zn = sqrt(b + 2*log(1/zn - 1/(zn*zn*zn) + 3/(zn*zn*zn*zn*zn)));
+      zn = sqrt(b + 2*log(((3/(zn*zn) - 1)/(zn*zn) + 1)/zn) );
+    }
+  
+    return zn;
+  }
+  
+  float Base2z::convertlogp2z(float logp) 
+    {
+      // logp must be the *natural* logarithm of p, not base 10
+      float z = 0.0;
+      
+      if(!issmalllogp(logp)) {
+	z = MISCMATHS::ndtri(exp(logp));
+      }
+      else {
+	z = logp2largez(logp);
+      }
+
+      return z;
+    }
+
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

base2z.h

+/*  base2z.h
+
+    Mark Woolrich & Mark Jenkinson, FMRIB Image Analysis Group
+
+    Copyright (C) 1999-2000 University of Oxford  */
+
+/*  CCOPYRIGHT  */
+
+#if !defined(__base2z_h)
+#define __base2z_h
+
+#include <iostream>
+#include <fstream>
+
+namespace MISCMATHS {
+  
+  class Base2z
+    {
+    public:
+      Base2z()
+	{}
+      
+      virtual ~Base2z() { delete base2z; }      
+      
+      float convertlogp2z(float logp); 
+      float logp2largez(float logp);
+      float logbeta(float v, float w);
+
+      virtual bool issmalllogp(float logp) = 0;
+
+    private:
+      
+      const Base2z& operator=(Base2z&);
+      Base2z(Base2z&);
+      
+      static Base2z* base2z;
+
+    };
+
+}
+
+#endif
+

f2z.cc

+/*  f2z.cc
+
+    Mark Woolrich & Mark Jenkinson, FMRIB Image Analysis Group
+
+    Copyright (C) 1999-2000 University of Oxford  */
+
+/*  CCOPYRIGHT  */
+
+#include <cmath>
+#include "f2z.h"
+#include "Log.h"
+#include "tracer_plus.h"
+#include <stdexcept>
+#include "libprob.h"
+
+using namespace NEWMAT;
+using namespace Utilities;
+
+namespace MISCMATHS {
+
+  F2z* F2z::f2z = NULL;
+ 
+  float F2z::largef2logp(float f, int d1, int d2)
+  {
+    Tracer_Plus ts("F2z::largef2logp");
+
+    // no of iterations:
+    int N = 20;
+
+    if (f<=0.0) {
+      cerr << "f cannot be zero or negative!" << endl;
+      return 0.0;
+    }
+  
+    if (d1<=0 || d2<=0) { 
+      cerr << "DOFs cannot be zero or negative!" << endl;
+      return 0.0; 
+    }
+ 
+    double alpha=d1/(double)d2;
+    double m=(d1+d2)/2.0;
+    double n=(1-d1/2.0);
+    double loggam = (d1/2.0)*(::log(d1/(double)d2)-logbeta(d2/2.0,d1/2.0));
+
+    //iter=f^(-n)/(alpha*(n+m-1)) + n*f^(-(n+1))/(alpha^2*(n+m-1)*(n+m)) + n*(n+1)*f^(-(n+2))/(alpha^3*(n+m-1)*(n+m)*(n+m+1));
+
+    double top = 1.0;
+    double bot = n+m-1;
+    double iter = 0.0;
+ 
+//     cerr << "logbeta(d2/2.0,d1/2.0)=" << logbeta(d2/2.0,d1/2.0) << endl;
+//     cerr << "loggam = " << loggam << endl;
+//     cerr << "n = " << n << endl;
+//     cerr << "m = " << m << endl;
+
+    for(int i = 1; i <= N; i++)
+      {
+	//	cerr << "i=" << i;
+	iter = iter + top*::pow(f,(-(n+i-1)))/(::pow(alpha,i)*bot);	
+	top = top*(n-1+i)*(-1);
+	bot = bot*(n+m-1+i);
+      }
+
+    float logp = loggam-(m-1)*(::log(1+alpha*f))+::log(iter);
+
+    //      cerr << "iter = " << iter << endl;
+    //      cerr << "logp = " << logp << endl;
+
+    return logp;
+  }
+
+  bool F2z::islargef(float f, int d1, int d2, float &logp) {
+   
+    if(f > 1.0 && d1>1)
+      {
+	logp=largef2logp(f,d1,d2);	
+	return issmalllogp(logp);
+      }
+    else
+      return false;
+  }
+
+  bool F2z::issmalllogp(float logp) {
+    return (logp < -14.5);
+  }
+
+  float F2z::convert(float f, int d1, int d2) 
+  {
+    Tracer_Plus ts("F2z::convert");
+
+    float z = 0.0, logp=0.0;
+
+    if(!islargef(f,d1,d2,logp)) {
+
+      double p = MISCMATHS::fdtr(d1, d2, f);
+
+      z = MISCMATHS::ndtri(p);
+    }
+      else {
+
+	z = logp2largez(logp);
+      }
+
+      return z;
+    }
+
+  void F2z::ComputeFStats(const ColumnVector& p_fs, int p_dof1, int p_dof2, ColumnVector& p_zs)
+  {
+    ColumnVector dof2 = p_fs;
+    dof2 = p_dof2;
+    ComputeFStats(p_fs,p_dof1,dof2,p_zs);
+  }
+  
+  void F2z::ComputeFStats(const ColumnVector& p_fs, int p_dof1, const ColumnVector& p_dof2, ColumnVector& p_zs)
+  {
+    Tracer_Plus ts("F2z::ComputeFStats");
+    
+    int numTS = p_fs.Nrows();
+    
+    p_zs.ReSize(numTS);
+    F2z& f2z = F2z::getInstance();
+    
+    for(int i = 1; i <= numTS; i++)
+      {
+	// 	  cerr << "i=" << i << endl;
+	// 	  cerr << "p_fs(i)=" << p_fs(i) << endl;
+	// 	  cerr << "p_dof1=" << p_dof1 << endl;
+	// 	  cerr << "p_dof2=" << p_dof2 << endl;
+	
+	if (p_fs(i) > 0.0)
+	  {
+	    p_zs(i) = f2z.convert(p_fs(i),p_dof1,p_dof2(i));  
+	  }
+	else
+	  {
+	    p_zs(i) = 0.0;
+	  }     
+      }
+  }
+  
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

f2z.h

+/*  f2z.h
+
+    Mark Woolrich & Mark Jenkinson, FMRIB Image Analysis Group
+
+    Copyright (C) 1999-2000 University of Oxford  */
+
+/*  CCOPYRIGHT  */
+
+#if !defined(__f2z_h)
+#define __f2z_h
+  
+#include <iostream>
+#include <fstream>
+#include "newmatap.h"
+#include "newmatio.h"
+#include "base2z.h"
+//#include "miscmaths.h"
+
+using namespace NEWMAT;
+
+namespace MISCMATHS {
+
+  class F2z : public Base2z
+    {
+    public:
+      static F2z& getInstance();
+      ~F2z() { delete f2z; }
+      
+      float convert(float f, int d1, int d2);
+
+      static void ComputeFStats(const ColumnVector& p_fs, int p_dof1, int p_dof2, ColumnVector& p_zs);
+      static void ComputeFStats(const ColumnVector& p_fs, int p_dof1, const ColumnVector& p_dof2, ColumnVector& p_zs);
+
+    private:
+      F2z() : Base2z()
+	{}
+      
+      const F2z& operator=(F2z&);
+      F2z(F2z&);
+      
+      bool issmalllogp(float logp);
+      bool islargef(float t, int d1, int d2, float &logp);
+      float largef2logp(float t, int d1, int d2);
+
+      static F2z* f2z;
+
+    };
+ 
+  inline F2z& F2z::getInstance(){
+    if(f2z == NULL)
+      f2z = new F2z();
+  
+    return *f2z;
+  }
+
+}
+
+#endif

histogram.cc

+/*  histogram.cc
+
+    Mark Woolrich, FMRIB Image Analysis Group
+
+    Copyright (C) 1999-2000 University of Oxford  */
+
+/*  CCOPYRIGHT  */
+
+#include "histogram.h"
+#include "miscmaths.h"
+
+#ifndef NO_NAMESPACE
+namespace MISCMATHS {
+#endif
+
+  void Histogram::generate()
+    {
+      Tracer ts("Histogram::generate");
+
+      int size = sourceData.Nrows();
+      
+      if(calcRange)
+	{
+	  // calculate range automatically
+	  histMin=histMax=sourceData(1);
+	  for(int i=1; i<=size; i++)
+	    {
+	      if (sourceData(i)>histMax)
+		histMax=sourceData(i);
+	      if (sourceData(i)<histMin)
+		histMin=sourceData(i);
+	    }
+	}
+      
+      // zero histogram
+      histogram.ReSize(bins);
+      histogram=0;
+      
+      // create histogram; the MIN is so that the maximum value falls in the
+      // last valid bin, not the (last+1) bin
+      for(int i=1; i<=size; i++)
+	{
+	   histogram(getBin(sourceData(i)))++;
+	}
+    }
+
+  int Histogram::integrate(float value1, float value2) const
+    {
+      int upperLimit = getBin(value2);
+      int sum = 0;
+
+      for(int i = getBin(value1)+1; i< upperLimit; i++)
+	{
+	  sum += (int)histogram(i);
+	}
+      return sum;
+    }
+
+  float Histogram::mode() const
+    {
+      int maxbin = 0;
+      int maxnum = 0;
+
+      for(int i = 1; i< bins; i++)
+	{
+	  if((int)histogram(i) > maxnum) {
+	    maxnum = (int)histogram(i);
+	    maxbin = i;
+	  }
+	}
+
+      return getValue(maxbin);
+    }
+
+#ifndef NO_NAMESPACE
+}
+#endif
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

histogram.h

+/*  histogram.h
+
+    Mark Woolrich, FMRIB Image Analysis Group
+
+    Copyright (C) 1999-2000 University of Oxford  */
+
+/*  CCOPYRIGHT  */
+
+#if !defined(__histogram_h)
+#define __histogram_h
+
+#include <iostream>
+#include <fstream>
+#define WANT_STREAM
+#define WANT_MATH
+
+#include "newmatap.h"
+#include "newmatio.h"
+#include "miscmaths.h"
+
+using namespace NEWMAT;
+
+namespace MISCMATHS {
+ 
+  class Histogram
+    {
+    public:
+      Histogram(const ColumnVector& psourceData, int numBins)
+	: sourceData(psourceData), calcRange(true), bins(numBins){}
+
+      Histogram(const ColumnVector& psourceData, float phistMin, float phistMax, int numBins) 
+	: sourceData(psourceData), calcRange(false), histMin(phistMin), histMax(phistMax), bins(numBins){}
+      
+      void generate();
+      
+      float getHistMin() const {return histMin;}
+      float getHistMax() const {return histMax;}
+      void setHistMax(float phistMax) {histMax = phistMax;}
+      void setHistMin(float phistMin) {histMin = phistMin;}
+
+      int integrateAll() {return sourceData.Nrows();}
+
+      const ColumnVector& getData() {return histogram;}
+
+      int integrateToInf(float value) const { return integrate(value, histMax); }
+      int integrateFromInf(float value) const { return integrate(histMin, value); }
+      int integrate(float value1, float value2) const;
+
+      float mode() const;
+
+      int getBin(float value) const;
+      float getValue(int bin) const;
+
+    protected:
+
+    private:
+      Histogram();
+      const Histogram& operator=(Histogram&);
+      Histogram(Histogram&);
+
+      const ColumnVector& sourceData;
+      ColumnVector histogram;
+
+      bool calcRange;
+
+      float histMin;
+      float histMax;
+
+      int bins; // number of bins in histogram
+    };
+
+  inline int Histogram::getBin(float value) const
+    {
+      return Max(1, Min((int)((((float)bins)*((float)(value-histMin)))/((float)(histMax-histMin))),bins));
+    }
+  
+  inline float Histogram::getValue(int bin) const
+    {
+      return (bin*(histMax-histMin))/(float)bins + histMin;
+    }
+
+}
+
+#endif
+
+
+
+
+

kernel.cc

+/*  kernel.cc
+
+    Mark Jenkinson, FMRIB Image Analysis Group
+
+    Copyright (C) 2001 University of Oxford  */
+
+/*  CCOPYRIGHT  */
+
+#include "kernel.h"
+#include "miscmaths.h"
+
+namespace MISCMATHS {
+
+  set<kernelstorage*, kernelstorage::comparer> kernel::existingkernels;  
+
+  //////// Support function /////////
+
+  float kernelval(float x, int w, const ColumnVector& kernel)
+  {
+    // linearly interpolates to get the kernel at the point (x)
+    //   given the half-width w
+    if (fabs(x)>w) return 0.0;
+    float halfnk = (kernel.Nrows()-1.0)/2.0;
+    float dn = x/w*halfnk + halfnk + 1.0;
+    int n = (int) floor(dn);
+    dn -= n;
+    if (n>(kernel.Nrows()-1)) return 0.0;
+    if (n<1) return 0.0;
+    
+    return kernel(n)*(1.0-dn) + kernel(n+1)*dn;
+  }
+  
+  inline bool in_bounds(ColumnVector data, int index)
+  { return ( (index>=1) && (index<=data.Nrows())); }
+  
+  inline bool in_bounds(ColumnVector data, float index)
+  { return ( ((int)floor(index)>=1) && ((int)ceil(index)<=data.Nrows())); }
+  
+  // Support Functions
+
+  float sincfn(float x)
+  {
+    if (fabs(x)<1e-7) { return 1.0-fabs(x); }
+    float y=M_PI*x;
+    return sin(y)/y;
+  }
+  
+  float hanning(float x, int w)
+  {  // w is half-width
+    if (fabs(x)>w) 
+      return 0.0;
+    else
+      return (0.5 + 0.5 *cos(M_PI*x/w));
+  }
+  
+  float blackman(float x, int w)
+  {  // w is half-width
+    if (fabs(x)>w) 
+      return 0.0;
+    else
+      return (0.42 + 0.5 *cos(M_PI*x/w) + 0.08*cos(2.0*M_PI*x/w));
+  }
+  
+  float rectangular(float x, int w)
+  {  // w is half-width
+    if (fabs(x)>w) 
+      return 0.0;
+    else
+      return 1.0;
+  }
+
+  ColumnVector sinckernel1D(const string& sincwindowtype, int w, int n)
+  {  // w is full-width
+    int nstore = n;
+    if (nstore<1) nstore=1;
+    ColumnVector ker(nstore);
+    int hw = (w-1)/2; // convert to half-width
+    // set x between +/- width
+    float halfnk = (nstore-1.0)/2.0;
+    for (int n=1; n<=nstore; n++) {
+      float x=(n-halfnk-1)/halfnk*hw;
+      if ( (sincwindowtype=="hanning") || (sincwindowtype=="h") ) {
+	ker(n) = sincfn(x)*hanning(x,hw);
+      } else if ( (sincwindowtype=="blackman") || (sincwindowtype=="b") ) {
+	ker(n) = sincfn(x)*blackman(x,hw);
+      } else if ( (sincwindowtype=="rectangular") || (sincwindowtype=="r") ) {
+	ker(n) = sincfn(x)*rectangular(x,hw);
+      } else {
+	cerr << "ERROR: Unrecognised sinc window type - using Blackman" << endl;
+	ker = sinckernel1D("b",w,nstore);
+	return ker;
+      }
+    }
+    return ker;
+  }
+
+  
+  kernel sinckernel(const string& sincwindowtype, int w, int nstore) 
+  {
+    kernel sinck;
+    sinck = sinckernel(sincwindowtype,w,w,w,nstore);
+    return sinck;
+  }
+
+
+  kernel sinckernel(const string& sincwindowtype,
+		    int wx, int wy, int wz, int nstore)
+  { // widths are full-widths
+    kernel sinckern;
+    if (nstore<1) nstore=1;
+
+    // convert all widths to half-widths
+    int hwx = (wx-1)/2;
+    int hwy = (wy-1)/2;
+    int hwz = (wz-1)/2;
+
+    ColumnVector kx, ky, kz;
+    // calculate kernels
+    kx = sinckernel1D(sincwindowtype,wx,nstore);
+    ky = sinckernel1D(sincwindowtype,wy,nstore);
+    kz = sinckernel1D(sincwindowtype,wz,nstore);
+    
+    sinckern.setkernel(kx,ky,kz,hwx,hwy,hwz);
+    return sinckern;
+  }
+
+  // dummy fn for now
+  float extrapolate_1d(const ColumnVector data, const int index)
+  {
+    float extrapval;
+
+    if (in_bounds(data, index))
+      extrapval = data(index);
+    else if (in_bounds(data, index-1))
+      extrapval = data(data.Nrows());
+    else if (in_bounds(data, index+1))
+      extrapval = data(1);
+    else
+      extrapval = mean(data).AsScalar();
+
+    return extrapval;
+  }
+  
+  // basic trilinear call
+  float interpolate_1d(ColumnVector data, const float index)
+  {
+    float interpval;
+    int low_bound = (int)floor(index);
+    int high_bound = (int)ceil(index); 
+
+    if (in_bounds(data, index))
+      interpval = data(low_bound) + (index - low_bound)*(data(high_bound) - data(low_bound));
+    else
+      interpval = extrapolate_1d(data, round(index));
+
+    return interpval;
+  }
+
+    
+  //////// Spline Support /////////
+
+  float hermiteinterpolation_1d(ColumnVector data, int p1, int p4, float t)
+  {
+    // Q(t) = (2t^3 - 3t^2 + 1)P_1 + (-2t^3 + 3t^2)P_4 + (t^3 - 2t^2 + t)R_1 + (t^3 - t^2)R_4
+    // inputs: points P_1, P_4; tangents R_1, R_4; interpolation index t;
+
+    float retval, r1, r4;
+    
+    if (!in_bounds(data,p1) || !in_bounds(data,p4)) {
+      cerr << "ERROR: One or more indicies lie outside the data range. Returning ZERO" << endl;
+      retval = 0.0;
+    } else if ((t < 0) || (t > 1)) {
+      cerr << "ERROR: Interpolation index must lie between 0 and 1. Returning ZERO" << endl;
+      retval = 0.0;
+      /*    } else if (t == 0.0) {
+	    retval = data(p1);
+	    } else if (t == 1.0) {
+	    retval = data(p4);
+      */   
+    } else {
+      r1 = 0.5 * (extrapolate_1d(data, p1) - extrapolate_1d(data, p1 - 1)) + 0.5 * (extrapolate_1d(data, p1 + 1) - extrapolate_1d(data, p1));// tangent @ P_1
+      r4 = 0.5 * (extrapolate_1d(data, p4) - extrapolate_1d(data, p4 - 1)) + 0.5 * (extrapolate_1d(data, p4 + 1) - extrapolate_1d(data, p4));// tangent @ P_4
+      
+      float t2 = t*t; float t3 = t2*t;
+      retval = (2*t3 - 3*t2 + 1)*data(p1) + (-2*t3 + 3*t2)*data(p4) + (t3 - 2*t2 + t)*r1 + (t3 - t2)*r4;
+    }
+
+    cerr << "p1, p4, t, r1, r4 = " << p1 << ", " << p4 << ", " << t << ", " << r1 << ", " << r4 << endl;
+    
+    return retval;
+  }
+
+
+  //////// Kernel Interpolation Call /////////
+
+  float kernelinterpolation_1d(ColumnVector data, float index, ColumnVector userkernel, int width)
+  {
+    int widthx = (width - 1)/2;
+    // kernel half-width  (i.e. range is +/- w)
+
+    int wx= widthx;
+    ColumnVector kernelx = userkernel;
+    float *storex = new float[2*wx+1];
+  
+    int ix0;
+    ix0 = (int) floor(index);
+    
+    float convsum=0.0, interpval=0.0, kersum=0.0;
+    
+    for (int d=-wx; d<=wx; d++) {
+      storex[d+wx] = kernelval((index-ix0+d),wx,kernelx);
+    }
+
+    int xj;
+    for (int x1=ix0-wx; x1<=ix0+wx; x1++) {
+      if (in_bounds(data, x1)) {
+	xj=ix0-x1+wx;
+	float kerfac = storex[xj];
+	//	cerr << "x1 = " << x1 << endl;
+	//      cerr << "data(x1) = " << data(x1) << endl;
+	convsum += data(x1) * kerfac;
+	kersum += kerfac;
+      }
+    }
+
+    delete(storex);
+    
+    if ( (fabs(kersum)>1e-9) ) {
+      interpval = convsum / kersum;
+    } else {
+      interpval = (float) extrapolate_1d(data, ix0);
+    }
+
+    // cerr << "interpval = " << interpval << endl;
+
+    return interpval;
+  }
+  
+  ////// Kernel wrappers //////
+
+  float kernelinterpolation_1d(ColumnVector data, float index)
+  {
+    ColumnVector userkernel = sinckernel1D("hanning", 7, 1201);
+    return kernelinterpolation_1d(data, index, userkernel, 7);
+  }
+
+  float kernelinterpolation_1d(RowVector data, float index)
+  {
+    ColumnVector userkernel = sinckernel1D("hanning", 7, 1201);
+    return kernelinterpolation_1d(data.t(), index, userkernel, 7);
+  }
+}
+
+
+
+
+
+
+
+
+
+
+
+ 

kernel.h

+/*  kernel.h
+
+    Mark Jenkinson, FMRIB Image Analysis Group
+
+    Copyright (C) 2001 University of Oxford  */
+
+/*  CCOPYRIGHT  */
+
+// General kernel interpolation class
+
+#if !defined(kernel_h)
+#define kernel_h
+
+#include <iostream>
+#include <string>
+#include <set>
+#include <cmath>
+#include "newmat.h"
+      
+using namespace NEWMAT;
+using namespace std;
+
+namespace MISCMATHS {
+
+  /////////////////////////////////////////////////////////////////////////
+
+  // Interpolation kernel storage class
+
+  class kernelstorage
+    {
+    private:
+      // NB: all widths are kernel half-widths (i.e. x \in [ -w, +w ] )
+      int p_widthx;
+      int p_widthy;
+      int p_widthz;
+      ColumnVector p_kernelx;
+      ColumnVector p_kernely;
+      ColumnVector p_kernelz;
+
+      // stop all forms of creation except the constructors below
+      kernelstorage();
+      const kernelstorage& operator=(kernelstorage&);
+      kernelstorage(kernelstorage&);
+
+
+    public:
+      float *storex;
+      float *storey;
+      float *storez;
+
+      kernelstorage(const ColumnVector& kx, const ColumnVector& ky,
+		    const ColumnVector& kz, int wx, int wy, int wz)	
+        {
+	  p_kernelx = kx; p_kernely = ky; p_kernelz = kz;
+	  p_widthx = wx;  p_widthy = wy;  p_widthz = wz;
+	  storez = new float[2*wz+1];
+	  storey = new float[2*wy+1];
+	  storex = new float[2*wx+1];
+ 	}
+
+      ~kernelstorage()
+      { 
+	delete storex;
+	delete storey;
+	delete storez;
+      }
+      
+      class comparer
+	{
+	public:
+	  bool operator()(const kernelstorage* k1, 
+			  const kernelstorage* k2) const
+	    {
+	      // comparison of sizes and values (toleranced)
+	      if ( (k1->p_widthx!=k2->p_widthx) || 
+		   (k1->p_widthy!=k2->p_widthy) || 
+		   (k1->p_widthz!=k2->p_widthz) )
+		return false;
+	      if ( ( (k1->p_kernelx - k2->p_kernelx).MaximumAbsoluteValue()
+		     > 1e-8 * k1->p_kernelx.MaximumAbsoluteValue() ) ||
+		   ( (k1->p_kernely - k2->p_kernely).MaximumAbsoluteValue() 
+		     > 1e-8 * k1->p_kernely.MaximumAbsoluteValue() ) ||
+		   ( (k1->p_kernelz - k2->p_kernelz).MaximumAbsoluteValue() 
+		     > 1e-8 * k1->p_kernelz.MaximumAbsoluteValue() ) )
+		return false;
+	      return true;
+	    }
+	};
+
+      friend class comparer;
+
+      int widthx() const { return p_widthx; }
+      int widthy() const { return p_widthy; }
+      int widthz() const { return p_widthz; }
+      const ColumnVector& kernelx() const { return p_kernelx; }
+      const ColumnVector& kernely() const { return p_kernely; }
+      const ColumnVector& kernelz() const { return p_kernelz; }
+
+    };
+
+
+  /////////////////////////////////////////////////////////////////////////////
+
+  class kernel
+    {
+    private:
+      static set<kernelstorage*, kernelstorage::comparer> existingkernels;
+      kernelstorage* storedkernel;
+
+    public:
+      kernel() { storedkernel = 0; }
+
+      const kernel& operator=(const kernel& source)
+      {
+	// am allowed to copy pointers if other class either
+	//  always exists or manages reference counts and self-deletes
+	this->existingkernels = source.existingkernels;
+	this->storedkernel = source.storedkernel;
+	// signal storedkernel has an extra reference
+	//   and that old storedkernel has one less reference
+	return *this;
+      }
+
+      kernel(const kernel& source)
+      {
+	this->operator=(source);
+      }
+
+      virtual ~kernel() 
+      { 
+	// signal storedkernel it has one less reference
+      }
+      
+      
+      void setkernel(const ColumnVector& kx, const ColumnVector& ky,
+		     const ColumnVector& kz, int wx, int wy, int wz)
+      {		  
+	// see if already in list:
+	storedkernel = new kernelstorage(kx,ky,kz,wx,wy,wz);
+	set<kernelstorage*, kernelstorage::comparer>::iterator 
+	  it = existingkernels.find(storedkernel);
+	if (it==existingkernels.end()) {		  
+	  existingkernels.insert(storedkernel);
+	  // signal that this is the first reference for storedkernel
+	} else {
+	  delete storedkernel;
+	  storedkernel = *it;
+	  // signal that *it has another reference now
+	}
+      }
+
+      const kernelstorage* kernelvals() { return storedkernel; }
+      
+  };
+
+
+  /////////////////////////////////////////////////////////////////////////
+
+  //////// Support functions /////////
+  
+  float kernelval(float x, int w, const ColumnVector& kernel);
+  float sincfn(float x);
+  float hanning(float x, int w);
+  float blackman(float x, int w);
+  float rectangular(float x, int w);
+  ColumnVector sinckernel1D(const string& sincwindowtype, int w, int n);
+  kernel sinckernel(const string& sincwindowtype, int w, int nstore);
+  kernel sinckernel(const string& sincwindowtype, 
+		    int wx, int wy, int wz, int nstore);
+  float extrapolate_1d(const ColumnVector data, const int index);
+  float interpolate_1d(ColumnVector data, const float index);
+  float kernelinterpolation_1d(ColumnVector data, float index, ColumnVector userkernel, int width);
+  float kernelinterpolation_1d(ColumnVector data, float index);
+  float kernelinterpolation_1d(RowVector data, float index);
+  float hermiteinterpolation_1d(ColumnVector data, int p1, int p4, float t);
+}
+
+#endif
+

kertest.cc

+#include "kernel.h"
+
+using namespace MISCMATHS;
+using namespace NEWMAT;
+
+int main (int argc,char** argv)
+{
+  cerr << "Test program for 1D kernel interpolation" << endl;
+
+  ColumnVector data(10);
+  data << 0 << 0 << 0 << 0 << 0 << 1 << 1 << 1 << 1 << 1;
+  ColumnVector newVec = data;
+  
+  cerr << "Input: " << data << endl;
+
+  for (int index = 0; index <= 10; index++)
+    newVec[index] = kernelinterpolation_1d(data, index+0.5);
+
+  cerr << "Result: " << newVec << endl;
+  
+  return 0;
+}

miscmaths.h

+/*  miscmaths.h
+
+    Mark Jenkinson & Mark Woolrich & Christian Beckmann, FMRIB Image Analysis Group
+
+    Copyright (C) 1999-2000 University of Oxford  */
+
+/*  CCOPYRIGHT  */
+
+// Miscellaneous maths functions
+
+
+#if !defined(__miscmaths_h)
+#define __miscmaths_h
+
+#include <cmath>
+#include <vector>
+#include <iostream>
+#include <assert.h>
+#include <fstream>
+#include <strstream>
+#include <string>
+#include "newmatap.h"
+#include "kernel.h"
+
+//#pragma interface
+
+using namespace NEWMAT;
+using namespace std;
+
+namespace MISCMATHS {
+
+#ifndef M_PI
+#define M_PI            3.14159265358979323846
+#endif
+
+#define OUT(t) cout<<#t "="<<t<<endl;
+
+  // IO/string stuff
+  template <class T> string num2str(T n, int width=-1);
+
+  bool isnum(const string& str);
+  ReturnMatrix read_ascii_matrix(const string& filename, int nrows, int ncols);
+  ReturnMatrix read_ascii_matrix(int nrows, int ncols, const string& filename);
+  ReturnMatrix read_ascii_matrix(const string& filename);
+  ReturnMatrix read_vest(string p_fname);
+  ReturnMatrix read_binary_matrix(const string& filename);
+
+  int write_ascii_matrix(const Matrix& mat, const string& filename, 
+			 int precision=-1);
+  int write_ascii_matrix(const string& filename, const Matrix& mat, 
+			 int precision=-1);
+  int write_vest(const Matrix& x, string p_fname, int precision=-1);
+  int write_vest(string p_fname, const Matrix& x, int precision=-1);
+  int write_binary_matrix(const Matrix& mat, const string& filename);
+
+  // more basic IO calls
+  string skip_alpha(ifstream& fs);
+  ReturnMatrix read_ascii_matrix(ifstream& fs, int nrows, int ncols);
+  ReturnMatrix read_ascii_matrix(int nrows, int ncols, ifstream& fs);
+  ReturnMatrix read_ascii_matrix(ifstream& fs);
+  ReturnMatrix read_binary_matrix(ifstream& fs);
+  int write_ascii_matrix(const Matrix& mat, ofstream& fs, int precision=-1);
+  int write_ascii_matrix(ofstream& fs, const Matrix& mat, int precision=-1);
+  int write_binary_matrix(const Matrix& mat, ofstream& fs);
+
+
+  // General maths
+
+  int round(int x); 
+  int round(float x);
+  int round(double x);
+
+  int sign(int x);
+  int sign(float x);
+  int sign(double x);
+
+  inline double pow(double x, float y) { return std::pow(x,(double) y); }
+  inline double pow(float x, double y) { return std::pow((double) x,y); }
+  inline double pow(double x, int y) { return std::pow(x,(double) y); }
+  inline float pow(float x, int y) { return std::pow(x,(float) y); }
+  inline double pow(int x, double y) { return std::pow((double)x, y); }
+  inline float pow(int x, float y) { return std::pow((float)x, y); }
+
+  inline double sqrt(int x) { return std::sqrt((double) x); }
+  inline double log(int x) { return std::log((double) x); }
+
+  int periodicclamp(int x, int x1, int x2);
+
+  template<class S, class T> 
+   inline T Min(const S &a, const T &b) { if (a<b) return (T) a; else return b; }
+
+  template<class S, class T>
+   inline T Max(const S &a, const T &b) { if (a>b) return (T) a; else return b; }
+
+  template<class T>
+   inline T Sqr(const T& x) { return x*x; } 
+
+  ColumnVector cross(const ColumnVector& a, const ColumnVector& b);
+  ColumnVector cross(const Real *a, const Real *b);
+
+  inline float dot(const ColumnVector& a, const ColumnVector& b)
+    { return Sum(SP(a,b)); }
+
+  double norm2(const ColumnVector& x);
+
+  int Identity(Matrix& m);
+  ReturnMatrix Identity(int num);
+
+  int diag(Matrix& m, const float diagvals[]);
+  int diag(Matrix& m, const ColumnVector& diagvals);
+  ReturnMatrix diag(const Matrix& mat);
+  ReturnMatrix pinv(const Matrix& mat);
+  ReturnMatrix sqrtaff(const Matrix& mat);
+
+  void reshape(Matrix& r, const Matrix& m, int nrows, int ncols);
+  ReturnMatrix reshape(const Matrix& m, int nrows, int ncols);
+
+  int construct_rotmat_euler(const ColumnVector& params, int n, Matrix& aff);
+  int construct_rotmat_euler(const ColumnVector& params, int n, Matrix& aff,
+		   const ColumnVector& centre);
+  int construct_rotmat_quat(const ColumnVector& params, int n, Matrix& aff);
+  int construct_rotmat_quat(const ColumnVector& params, int n, Matrix& aff,
+		   const ColumnVector& centre);
+  int make_rot(const ColumnVector& angl, const ColumnVector& centre, 
+	       Matrix& rot);
+
+  int getrotaxis(ColumnVector& axis, const Matrix& rotmat);
+  int rotmat2euler(ColumnVector& angles, const Matrix& rotmat);
+  int rotmat2quat(ColumnVector& quaternion, const Matrix& rotmat);
+  int decompose_aff(ColumnVector& params, const Matrix& affmat,
+		    int (*rotmat2params)(ColumnVector& , const Matrix& ));
+  int decompose_aff(ColumnVector& params, const Matrix& affmat,
+		    const ColumnVector& centre,
+		    int (*rotmat2params)(ColumnVector& , const Matrix& ));
+  int compose_aff(const ColumnVector& params, int n, const ColumnVector& centre,
+		  Matrix& aff, 
+		  int (*params2rotmat)(const ColumnVector& , int , Matrix& ,
+				       const ColumnVector& ) );
+  float rms_deviation(const Matrix& affmat1, const Matrix& affmat2, 
+		      const ColumnVector& centre, const float rmax); 
+  float rms_deviation(const Matrix& affmat1, const Matrix& affmat2, 
+		      const float rmax=80.0); 
+  float median(const ColumnVector& x);
+
+  // returns the first P such that 2^P >= abs(N). 
+  int nextpow2(int n);
+
+  // Auto-correlation function estimate of columns of p_ts
+  // gives unbiased estimate - scales the raw correlation by 1/(N-abs(lags))
+  void xcorr(const Matrix& p_ts, Matrix& ret, int lag = 0, int p_zeropad = 0);
+  ReturnMatrix xcorr(const Matrix& p_ts, int lag = 0, int p_zeropad = 0);
+
+  // removes trend from columns of p_ts
+  // if p_level==0 it just removes the mean
+  // if p_level==1 it removes linear trend
+  // if p_level==2 it removes quadratic trend
+  void detrend(Matrix& p_ts, int p_level=1);
+
+  ReturnMatrix zeros(const int dim1, const int dim2 = -1);
+  ReturnMatrix ones(const int dim1, const int dim2 = -1);
+  ReturnMatrix repmat(const Matrix& mat, const int rows = 1, const int cols = 1);
+  ReturnMatrix dist2(const Matrix& mat1, const Matrix& mat2);
+  ReturnMatrix abs(const Matrix& mat);
+  ReturnMatrix sqrt(const Matrix& mat);
+  ReturnMatrix log(const Matrix& mat);
+  ReturnMatrix exp(const Matrix& mat);
+  ReturnMatrix tanh(const Matrix& mat);
+  ReturnMatrix pow(const Matrix& mat, const double exp);
+  ReturnMatrix sum(const Matrix& mat, const int dim);
+  ReturnMatrix mean(const Matrix& mat, const int dim = 1);
+  ReturnMatrix var(const Matrix& mat, const int dim = 1);
+  ReturnMatrix max(const Matrix& mat);
+  ReturnMatrix max(const Matrix& mat,ColumnVector& index);
+  ReturnMatrix min(const Matrix& mat);
+  ReturnMatrix gt(const Matrix& mat1,const Matrix& mat2); 
+  ReturnMatrix lt(const Matrix& mat1,const Matrix& mat2); 
+  ReturnMatrix geqt(const Matrix& mat1,const Matrix& mat2); 
+  ReturnMatrix leqt(const Matrix& mat1,const Matrix& mat2); 
+  ReturnMatrix eq(const Matrix& mat1,const Matrix& mat2); 
+  ReturnMatrix neq(const Matrix& mat1,const Matrix& mat2); 
+
+
+   
+
+  void remmean(const Matrix& mat, Matrix& demeanedmat, Matrix& Mean,  const int dim = 1);
+  ReturnMatrix remmean(const Matrix& mat, const int dim = 1);
+  ReturnMatrix stdev(const Matrix& mat, const int dim = 1);
+  ReturnMatrix cov(const Matrix& mat, const int norm = 0);
+  ReturnMatrix corrcoef(const Matrix& mat, const int norm = 0);
+  void symm_orth(Matrix &Mat);
+  void powerspectrum(const Matrix &Mat1, Matrix &Result, bool useLog);
+
+  
+  ///////////////////////////////////////////////////////////////////////////
+  ///////////////////////////////////////////////////////////////////////////
+ 
+  // TEMPLATE DEFINITIONS //
+  
+  template <class T>
+  string num2str(T n, int width)
+  {
+    ostrstream os;
+    if (width>0) {
+      os.fill('0');
+      os.width(width);
+      os.setf(ios::internal, ios::adjustfield);
+    }
+    os << n << '\0';
+    return os.str();
+  }
+
+}
+
+#endif

miscprob.cc

+/*  miscprob.cc
+
+    Christian Beckmann & Mark Woolrich, FMRIB Image Analysis Group
+
+    Copyright (C) 1999-2000 University of Oxford  */
+
+/*  CCOPYRIGHT  */
+
+// Miscellaneous maths functions that rely on libprob
+
+#include "miscprob.h"
+#include "stdlib.h"
+#include "newmatio.h"
+#include <iostream>
+
+using namespace NEWMAT;
+
+namespace MISCMATHS {
+
+ReturnMatrix unifrnd(const int dim1, const int dim2, const float start, const float end)
+{
+  int tdim = dim2;
+  double tmpD=1.0;
+  if(tdim<0){tdim=dim1;}
+  Matrix res(dim1,tdim);
+  for (int mc=1; mc<=res.Ncols(); mc++) {
+    for (int mr=1; mr<=res.Nrows(); mr++) {
+      //tmpD = (rand()+1)/double(RAND_MAX+2.0);
+      //res(mr,mc)=(tmpD)*(end-start)+start;
+      drand(&tmpD);
+      res(mr,mc)=(tmpD-1)*(end-start)+start;
+    }
+  }
+
+  res.Release();
+  return res;
+}
+
+ReturnMatrix normrnd(const int dim1, const int dim2, const float mu, const float sigma)
+{
+  int tdim = dim2;
+  double tmpD=1.0;
+  if(tdim<0){tdim=dim1;}
+  Matrix res(dim1,tdim);
+  for (int mc=1; mc<=res.Ncols(); mc++) {
+    for (int mr=1; mr<=res.Nrows(); mr++) {
+      //tmpD = (rand()+1)/double(RAND_MAX+2.0);
+      //res(mr,mc)=ndtri(tmpD)*sigma+mu ;
+      drand(&tmpD);
+      res(mr,mc)=ndtri(tmpD-1)*sigma+mu ;
+    }
+  }
+
+  res.Release();
+
+  return res;
+}
+
+ReturnMatrix normpdf(const RowVector& vals, const float mu, const float sigma)
+{
+  RowVector res(vals);
+  for (int mc=1; mc<=res.Ncols(); mc++){
+    res(mc) = std::exp(-0.5*(std::pow(vals(mc)-mu,2)/sigma))*std::pow(2*M_PI*sigma,-0.5);
+  }
+
+  res.Release();
+  return res;
+}
+
+ReturnMatrix normcdf(const RowVector& vals, const float mu, const float var)
+{
+  RowVector res(vals);
+  RowVector tmp;
+  tmp = (vals-mu)/std::sqrt(var);
+  for (int mc=1; mc<=res.Ncols(); mc++){
+    res(mc) = ndtr(tmp(mc));
+  }
+
+  res.Release();
+  return res;
+}
+
+ReturnMatrix gammacdf(const RowVector& vals, const float mu, const float sigma)
+{
+  RowVector res(vals);
+  res=0;
+  if((mu>0)&&(sigma>0)){
+    float b = std::pow(mu,2)/sigma;
+    float a = mu/sigma;  
+    for (int mc=1; mc<=res.Ncols(); mc++){
+      if(vals(mc)>0)
+	res(mc) = gdtr(a,b,vals(mc));
+    }
+  }
+  res.Release();
+  return res;
+}
+
+ReturnMatrix gammapdf(const RowVector& vals, const float mu, const float sigma)
+{
+  RowVector res(vals);
+  res=0;
+  if((mu>0)&&(sigma>0.00001)){
+    float a = std::pow(mu,2)/sigma;
+    float b = mu/sigma;
+    float c = lgam(a);
+    if(std::abs(c) < 150){
+      for (int mc=1; mc<=res.Ncols(); mc++){
+	if(vals(mc)>0.000001){
+	  res(mc) = std::exp(a*std::log(b) + 
+			     (a-1) * std::log(vals(mc)) 
+			     - b*vals(mc) - c);
+	}
+      }
+    }
+  }
+  res.Release();
+  return res;
+}
+
+
+ReturnMatrix normpdf(const RowVector& vals, const RowVector& mu, const RowVector& sigma)
+{
+  Matrix res(mu.Ncols(),vals.Ncols());
+  for (int mc=1; mc<=res.Ncols(); mc++){
+    for (int mr=1; mr<=res.Nrows(); mr++){
+      res(mr,mc) = std::exp(-0.5*(std::pow(vals(mc)-mu(mr),2)/sigma(mr)))*std::pow(2*M_PI*sigma(mr),-0.5);
+    }
+  }
+
+  res.Release();
+  return res;
+}
+
+
+ReturnMatrix mvnormrandm(const RowVector& mu, const SymmetricMatrix& covar, int nsamp) 
+{     
+//   Matrix eig_vec; 
+//   DiagonalMatrix eig_val;
+//   EigenValues(covar,eig_val,eig_vec);
+
+//   Matrix ret = ones(nsamp, 1)*mu + dnormrandm(nsamp,mu.Ncols())*sqrt(eig_val)*eig_vec.t();
+  Mvnormrandm mvn(mu, covar);
+
+  return mvn.next(nsamp);
+}
+
+}

miscprob.h

+/*  miscprob.h
+
+    Christian Beckmann & Mark Woolrich, FMRIB Image Analysis Group
+
+    Copyright (C) 1999-2000 University of Oxford  */
+
+/*  CCOPYRIGHT  */
+
+// Miscellaneous maths functions that rely on libprob build ontop of miscmaths
+
+
+#if !defined(__miscprob_h)
+#define __miscprob_h
+
+#include "miscmaths.h"
+#include "libprob.h"
+#include "stdlib.h"
+
+using namespace NEWMAT;
+
+namespace MISCMATHS {
+ 
+  ReturnMatrix unifrnd(const int dim1 = 1, const int dim2 = -1, 
+		       const float start = 0, const float end = 1);
+  
+  ReturnMatrix normrnd(const int dim1 = 1, const int dim2 = -1, 
+		       const float mu = 0, const float sigma = 1);
+
+  ReturnMatrix mvnrnd(const RowVector& mu, const SymmetricMatrix& covar, int nsamp = 1);
+  
+
+  ReturnMatrix normpdf(const RowVector& vals, const float mu = 0, const float sigma = 1);
+
+  ReturnMatrix normpdf(const RowVector& vals, const RowVector& mus, 
+		       const RowVector& sigmas);
+
+  ReturnMatrix normcdf(const RowVector& vals, const float mu = 0, const float sigma = 1);
+
+  ReturnMatrix gammapdf(const RowVector& vals, const float mu = 0, const float sigma = 1);
+
+  ReturnMatrix gammacdf(const RowVector& vals, const float mu = 0, const float sigma = 1);
+
+  
+  class Mvnormrandm
+    {
+    public:
+      Mvnormrandm(){}
+
+      Mvnormrandm(const RowVector& pmu, const SymmetricMatrix& pcovar) :
+	mu(pmu),
+	covar(pcovar)
+	{
+	  Matrix eig_vec;
+	  DiagonalMatrix eig_val;
+	  EigenValues(covar,eig_val,eig_vec);
+
+	  covarw = sqrt(eig_val)*eig_vec.t();
+	}
+
+      ReturnMatrix next(int nsamp = 1) const 
+	{
+	  Matrix ret = ones(nsamp, 1)*mu + normrnd(nsamp,mu.Ncols())*covarw;
+	  ret.Release();
+	  return ret;
+	}
+    private:      
+
+      RowVector mu;
+      SymmetricMatrix covar;
+
+      Matrix covarw;
+
+    };
+}
+#endif
+
+
+
+
+
+

optimise.cc

+/*  optimise.cc
+
+    Mark Jenkinson, FMRIB Image Analysis Group
+
+    Copyright (C) 1999-2000 University of Oxford  */
+
+/*  CCOPYRIGHT  */
+
+// Mathematical optimisation functions
+
+
+#include <iostream>
+#include <cmath>
+#include "optimise.h"
+#include "miscmaths.h"
+
+namespace MISCMATHS {
+
+  // The following lines are ignored by the current SGI compiler
+  //  (version egcs-2.91.57)
+  // A temporary fix of including the std:: in front of all abs() etc
+  //  has been done for now
+  using std::abs;
+  using std::sqrt;
+  using std::exp;
+  using std::log;
+
+  bool estquadmin(float &xnew, float x1, float xmid, float x2, 
+		   float y1, float ymid, float y2)
+  {
+    // Finds the estimated quadratic minimum's position
+    float ad=0.0, bd=0.0, det=0.0;
+    ad = (xmid - x2)*(ymid - y1) - (xmid - x1)*(ymid - y2);
+    bd = -(xmid*xmid - x2*x2)*(ymid - y1) + (xmid*xmid - x1*x1)*(ymid - y2);
+    det = (xmid - x2)*(x2 -x1)*(x1 - xmid);
+    if ((fabs(det)>1e-15) && (ad/det < 0)) {  // quadratic only has a maxima
+      xnew = 0.0;
+      return false;
+    }
+    if (fabs(ad)>1e-15) {
+      xnew = -bd/(2*ad);
+      return true;
+    } else {  // near linear condition -> get closer to an end point
+      xnew = 0.0;
+      return false;
+    }
+    return false;
+  }
+
+
+  float extrapolatept(float x1, float xmid, float x2)
+  {
+    // xmid must be between x1 and x2
+    // use the golden ratio (scale similar result)
+    const float extensionratio = 0.3819660;
+    float xnew;
+    if (fabs(x2-xmid)>fabs(x1-xmid)) {
+      xnew = extensionratio * x2 + (1 - extensionratio) * xmid;
+    } else {
+      xnew = extensionratio * x1 + (1 - extensionratio) * xmid;
+    }
+    return xnew;
+  }
+  
+
+
+  float nextpt(float x1, float xmid, float x2, float y1, float ymid, float y2)
+  {
+    // x1 and x2 are the bounds, xmid is between them
+
+    float xnew;
+    bool quadok=false;
+    quadok = estquadmin(xnew,x1,xmid,x2,y1,ymid,y2);
+
+    // check to see that the quadratic result is in the range
+    if ((!quadok) || (xnew < Min(x1,x2)) || (xnew > Max(x1,x2))) {
+      xnew = extrapolatept(x1,xmid,x2);
+    }
+    return xnew;
+  }
+
+      
+
+  void findinitialbound(float &x1, float &xmid, float &x2, 
+			float &y1, float &ymid, float &y2, 
+			float (*func)(const ColumnVector &),
+			const ColumnVector &unitdir, const ColumnVector &pt)
+  {
+    const float extrapolationfactor = 1.6;
+    const float maxextrap = extrapolationfactor*2;
+    if (y1==0)  y1 = (*func)(x1*unitdir + pt);
+    if (ymid==0)  ymid = (*func)(xmid*unitdir + pt);
+    if (y1<ymid) {   // swap a and b if this is the case
+      float tempx = x1, tempy = y1;
+      x1 = xmid;     y1 = ymid;
+      xmid = tempx;  ymid = tempy;
+    }
+
+    float newx2 = 0.0, newy2=0.0, maxx2=0.0;
+    float dir=1.0;
+    if (xmid<x1) dir=-1.0;
+
+    bool quadok;
+
+    x2 = xmid + extrapolationfactor*(xmid - x1);
+    y2 = (*func)(x2*unitdir + pt);
+
+    while (ymid > y2) {  // note: must maintain y1 >= ymid
+	
+      // cout << "    <" << Min(x1,x2) << "," << xmid 
+      //   << "," << Max(x1,x2) << ">" << endl;
+      maxx2 = xmid + maxextrap*(x2 - xmid);
+      quadok = estquadmin(newx2,x1,xmid,x2,y1,ymid,y2);
+      if ((!quadok) || ((newx2 - x1)*dir<0) || ((newx2 - maxx2)*dir>0)) {
+	newx2 = xmid + extrapolationfactor*(x2-x1);
+      }
+      
+      newy2 = (*func)(newx2*unitdir + pt);
+
+      if ((newx2 - xmid)*(newx2 - x1)<0) {  // newx2 is between x1 and xmid
+	if (newy2 < ymid) {  // found a bracket!
+	  x2 = xmid;  y2 = ymid;
+	  xmid = newx2;  ymid = newy2;
+	  break;
+	} else {  // can use newx2 as a new value for x1 (as newy2 >= ymid)
+	  x1 = newx2;  y1 = newy2;
+	}
+      } else {  // newx2 is between xmid and maxx2
+	if (newy2 > ymid) { // found a bracket!
+	  x2 = newx2;  y2 = newy2;
+	  break;
+	} else if ((newx2 - x2)*dir<0) {  // newx2 closer to xmid than old x2
+	  x1 = xmid;  y1 = ymid;
+	  xmid = newx2;  ymid = newy2;
+	} else {
+	  x1 = xmid;  y1 = ymid;
+	  xmid = x2;  ymid = y2;
+	  x2 = newx2;  y2 = newy2;
+	}
+      }
+	
+    }
+
+    if ( (y2<ymid) || (y1<ymid) ) {
+      cerr << "findinitialbound failed to bracket: current triplet is" << endl;
+    }
+  }
+  
+
+  float optimise1d(ColumnVector &pt, const ColumnVector dir, 
+		  const ColumnVector tol, int &iterations_done, 
+		  float (*func)(const ColumnVector &), int max_iter,
+		  float init_value, float boundguess) 
+  {
+    // Golden Search Routine
+    // Must pass in the direction vector in N-space (dir), the initial
+    //  N-dim point (pt), the acceptable tolerance (tol) and other
+    //  stuff
+    // Note that the length of the direction vector is unimportant
+
+    float y1,y2,ymid;
+    float x1,x2,xmid;
+
+    // Calculate dot product of dir by tol
+    //  st (x1-x2)*dir_tol = average number of tolerances between x1 and x2
+    float dir_tol = 0.0;
+    ColumnVector unitdir;
+    unitdir = dir/std::sqrt(dir.SumSquare());
+    for (int n=1; n<=tol.Nrows(); n++) {
+      if (fabs(tol(n))>1e-15) {
+	dir_tol += fabs(unitdir(n)/tol(n));
+      }
+    }
+    float unittol = fabs(1/dir_tol);
+
+    // set up initial points
+    xmid = 0.0;
+    x1 = boundguess * unittol;  // initial guess (bound)
+    if (init_value==0.0) ymid = (*func)(xmid*unitdir + pt);
+    else ymid = init_value;
+    y1 = (*func)(x1*unitdir + pt);
+    findinitialbound(x1,xmid,x2,y1,ymid,y2,func,unitdir,pt);
+
+    // cout << "(" << x1 << "," << y1 << ")  ";
+    // cout << "(" << xmid << "," << ymid << ")  ";
+    // cout << "(" << x2 << "," << y2 << ")" << endl;
+
+    float min_dist = 0.1 * unittol;
+    float xnew, ynew;
+    int it=0;
+    while ( ((++it)<=max_iter) && (fabs((x2-x1)/unittol)>1.0) )
+      {
+	// cout << "  [" << Min(x1,x2) << "," << Max(x1,x2) << "]" << endl;
+
+	if (it>0) {
+	  xnew = nextpt(x1,xmid,x2,y1,ymid,y2);
+	} else {
+	  xnew = extrapolatept(x1,xmid,x2);
+	}
+
+	float dirn=1.0;
+	if (x2<x1) dirn=-1.0;
+
+	if (fabs(xnew - x1)<min_dist) {
+	  xnew = x1 + dirn*min_dist;
+	}
+
+	if (fabs(xnew - x2)<min_dist) {
+	  xnew = x2 - dirn*min_dist;
+	}
+
+	if (fabs(xnew - xmid)<min_dist) {
+	  xnew = extrapolatept(x1,xmid,x2);
+	}
+
+	if (fabs(xmid - x1)<0.4*unittol) {
+	  xnew = xmid + dirn*0.5*unittol;
+	}
+
+	if (fabs(xmid - x2)<0.4*unittol) {
+	  xnew = xmid - dirn*0.5*unittol;
+	}
+
+	ynew = (*func)(xnew*unitdir + pt);
+
+	if ((xnew - xmid)*(x2 - xmid) > 0) {  // is xnew between x2 and xmid ?
+	  // swap x1 and x2 so that xnew is between x1 and xmid
+	  float xtemp = x1;  x1 = x2;  x2 = xtemp;
+	  float ytemp = y1;  y1 = y2;  y2 = ytemp;
+	}
+	if (ynew < ymid) {
+	  // new interval is [xmid,x1] with xnew as best point in the middle
+	  x2 = xmid;  y2 = ymid;
+	  xmid = xnew;  ymid = ynew;
+	} else {
+	  // new interval is  [x2,xnew] with xmid as best point still
+	  x1 = xnew;  y1 = ynew;
+	}
+      }
+    iterations_done = it;
+    pt = xmid*unitdir + pt;
+    return ymid;
+  }
+
+
+
+  
+  float optimise(ColumnVector &pt, int numopt, const ColumnVector &tol, 
+		 float (*func)(const ColumnVector &), int &iterations_done, 
+		 int max_iter, const ColumnVector& boundguess)
+  {
+    // Calculate dot product of dir by tol
+    //  st (x1-x2)*dir_tol = average number of tolerances between x1 and x2
+    ColumnVector inv_tol(tol.Nrows());
+    inv_tol = 0.0;
+    for (int n=1; n<=tol.Nrows(); n++) {
+      if (fabs(tol(n))>1e-15) {
+	inv_tol(n) = fabs(1.0/tol(n));
+      }
+    }
+    inv_tol /= (float) tol.Nrows();
+
+    ColumnVector dir(pt.Nrows()), initpt;
+    int lit=0, littot=0, it=0;
+    float fval=0.0, bndguess;
+    while ((++it)<=max_iter)
+      {
+	initpt = pt;
+	bndguess = boundguess(Min(it,boundguess.Nrows()));  // ceiling of nrows
+	for (int n=1; n<=numopt; n++) {
+	  dir = 0.0;
+	  dir(n) = 1.0;
+	  fval = optimise1d(pt,dir,tol,lit,func,100,fval,bndguess);
+	  littot += lit;
+	}
+	// check to see if the point has moved more than one average tolerance
+	float avtol = SP((initpt - pt),inv_tol).SumAbsoluteValue();
+	if (avtol < 1.0) break;
+      }
+    // cout << endl << "Major iterations = " << it << endl;
+    iterations_done = littot;
+    return fval;
+  }
+
+
+}
+
+
+
+
+
+
+
+

optimise.h

+/*  optimise.h
+
+    Mark Jenkinson, FMRIB Image Analysis Group
+
+    Copyright (C) 1999-2000 University of Oxford  */
+
+/*  CCOPYRIGHT  */
+
+// Mathematical optimisation functions
+
+
+#if !defined(__optimise_h)
+#define __optimise_h
+
+#include <cmath>
+#include "newmatap.h"
+
+using namespace NEWMAT;
+
+namespace MISCMATHS {
+
+float optimise1d(ColumnVector &pt, const ColumnVector dir, 
+		const ColumnVector tol, int &iterations_done, 
+		float (*func)(const ColumnVector &), int max_iter,
+		float init_value, float boundguess);
+
+
+ float optimise(ColumnVector &pt, int numopt, const ColumnVector &tol, 
+		float (*func)(const ColumnVector &), int &iterations_done, 
+		int max_iter, const ColumnVector& boundguess);
+
+}
+
+#endif

t2z.cc

+/*  t2z.cc
+
+    Mark Woolrich & Mark Jenkinson, FMRIB Image Analysis Group
+
+    Copyright (C) 1999-2000 University of Oxford  */
+
+/*  CCOPYRIGHT  */
+
+#include <cmath>
+#include "t2z.h"
+#include "newmat.h"
+#include "tracer_plus.h"
+#include "libprob.h"
+
+using namespace NEWMAT;
+using namespace Utilities;
+
+namespace MISCMATHS {
+
+  T2z* T2z::t2z = NULL;
+ 
+  float T2z::larget2logp(float t, int dof)
+  {
+    //    static float logbeta[] = {   1.144729885849, 0.693147180560,
+    //			 0.451582705289, 0.287682072452, 
+    //			 0.163900632838, 0.064538521138,
+    //			 -0.018420923956, -0.089612158690, 
+    //			 -0.151952316581, -0.207395194346 } ;
+
+    //static const float pi = 3.141592653590;
+    //    static const float log2pi = log(2*pi);
+
+    // Large T extrapolation routine for converting T to Z values
+    //  written by Mark Jenkinson, March 2000
+    //
+    // It does T to Z via log(p) rather than p, since p becomes very
+    //  small and underflows the arithmetic
+    // Equations were derived by using integration by parts and give the
+    //  following formulae:
+    //  (1) T to log(p)    NB: n = DOF
+    //       log(p) = -1/2*log(n) - log(beta(n/2,1/2)) - (n-1)/2*log(1+t*t/n)
+    //                + log(1 - (n/(n+2))/(t*t) + 3*n*n/((n+2)*(n+4)*t*t*t*t)) 
+    //  (2) Z to log(p)
+    //       log(p) = -1/2*z*z - 1/2*log(2*pi) - log(z) 
+    //                + log(1 - 1/(z*z) + 3/(z*z*z*z))
+    // equation (2) is then solved by the recursion:
+    //   z_0 = sqrt(2*(-log(p) - 1/2*log(2*pi)))
+    //   z_{n+1} = sqrt(2*(-log(p) - 1/2*log(2*pi) - log(z_n) 
+    //             + log(1 - 1/(zn*zn) + 3/(zn*zn*zn*zn))
+    // In practice this recursion is quite accurate in 3 to 5 iterations
+    // Equation (1) is accurate to 1 part in 10^3 for T>7.5  (any n)
+    // Equation (2) is accurate to 1 part in 10^3 for Z>3.12  (3 iterations)
+
+
+    if (t<0) {
+      return larget2logp(-t,dof);
+    }
+  
+    float logp, lbeta;
+
+    if (dof<=0) { 
+      cerr << "DOF cannot be zero or negative!" << endl;
+      return 0.0; 
+    }
+
+    float n = (float) dof;
+
+    // complete Beta function
+    lbeta = this->logbeta(1/2.0,n/2.0);
+    //if (dof<=10) {
+    //lbeta = logbeta[dof-1];
+    //} else {
+    //lbeta = log2pi/2 - log(n)/2 + 1/(4*n);
+    //}
+    
+    // log p from t value
+    // logp = log( (1 - n/((n+2)*t*t) + 3*n*n/((n+2)*(n+4)*t*t*t*t))/(sqrt(n)*t))
+    //          - ((n-1)/2)*log(1 + t*t/n) - lbeta;
+    logp = log(( (3*n*n/((n+2)*(n+4)*t*t) - n/(n+2))/(t*t) + 1)/(sqrt(n)*t))
+      - ((n-1)/2)*log(1 + t*t/n) - lbeta;
+
+    return logp;
+  }
+
+  bool T2z::islarget(float t, int dof, float &logp) {
+    // aymptotic formalae are valid if 
+    //   log(p) < -14.5  (derived from Z-statistic approximation error)
+    // For dof>=15, can guarantee that log(p)>-33 (p > 1e-14) if T<7.5
+    //  and so in this region use conventional means, not asymptotic
+    if ((dof>=15) && (fabs(t)<7.5)) { return false; }
+    logp=larget2logp(t,dof);
+    if (dof>=15)  return true;  // force asymptotic calc for all T>=7.5, D>=15
+    return issmalllogp(logp);
+  }
+
+  bool T2z::issmalllogp(float logp) {
+	// aymptotic formula accurate to 1 in 10^3 for Z>4.9
+	// which corresponds to log(p)=-14.5
+	return (logp < -14.5);
+      }
+
+  float T2z::convert(float t, int dof) 
+    {
+      float z = 0.0, logp=0.0;
+      
+      if(!islarget(t,dof,logp)) {
+	//	cerr << "t = " << t << endl;
+	double p = MISCMATHS::stdtr(dof, t);
+	//cerr << "p = " << p << endl;
+	z = MISCMATHS::ndtri(p);
+      }
+      else {
+	z = logp2largez(logp);
+	if (t<0) z=-z;
+      }
+
+      return z;
+    }
+
+
+  float T2z::converttologp(float t, int dof) 
+    {
+      float logp=0.0;
+      
+      if(!islarget(t,dof,logp)) {
+	//	cerr << "t = " << t << endl;
+	logp = log(MISCMATHS::stdtr(dof, t));
+      }
+
+      return logp;
+    }
+
+  void T2z::ComputePs(const ColumnVector& p_vars, const ColumnVector& p_cbs, int p_dof, ColumnVector& p_ps)
+    {
+      Tracer ts("T2z::ComputePs");
+
+      int numTS = p_vars.Nrows();
+
+      T2z& t2z = T2z::getInstance();
+
+      p_ps.ReSize(numTS);
+
+      for(int i = 1; i <= numTS; i++)
+	{
+	  //cerr << "var = " << p_vars(i) << " at index "<< i << endl;
+	  //cerr << "cb = " << p_cbs(i) << " at index "<< i << endl;
+	  if ( (p_vars(i) != 0.0) && (p_cbs(i) != 0.0) )
+	    {
+	      if(p_vars(i) < 0.0)
+		{
+		  //cerr << "var = " << p_vars(i) << " at index "<< i << endl;
+		  p_ps(i) = 0.0;
+		}
+	      else
+		{
+		  p_ps(i) = t2z.converttologp(p_cbs(i)/sqrt(p_vars(i)),p_dof);
+		  
+		  //if(p_zs(i) == 0.0)
+		  //cerr << " at index " << i << endl;
+		}
+	    }
+	  else
+	    p_ps(i) = 0.0;
+	}
+    }
+
+  void T2z::ComputeZStats(const ColumnVector& p_vars, const ColumnVector& p_cbs, int p_dof, ColumnVector& p_zs)
+    {
+      ColumnVector dof = p_vars;
+      dof = p_dof;
+      ComputeZStats(p_vars,p_cbs,dof,p_zs);
+    }
+
+  void T2z::ComputeZStats(const ColumnVector& p_vars, const ColumnVector& p_cbs, const ColumnVector& p_dof, ColumnVector& p_zs)
+    {
+      Tracer ts("T2z::ComputeStats");
+
+      int numTS = p_vars.Nrows();
+
+      T2z& t2z = T2z::getInstance();
+
+      p_zs.ReSize(numTS);
+
+      for(int i = 1; i <= numTS; i++)
+	{
+	  //cerr << "var = " << p_vars(i) << " at index "<< i << endl;
+	  //cerr << "cb = " << p_cbs(i) << " at index "<< i << endl;
+	  if ( (p_vars(i) != 0.0) && (p_cbs(i) != 0.0) )
+	    {
+	      if(p_vars(i) < 0.0)
+		{
+		  //cerr << "var = " << p_vars(i) << " at index "<< i << endl;
+		  p_zs(i) = 0.0;
+		}
+	      else
+		{
+		  p_zs(i) = t2z.convert(p_cbs(i)/sqrt(p_vars(i)),p_dof(i));
+		  
+		  //if(p_zs(i) == 0.0)
+		  //cerr << " at index " << i << endl;
+		}
+	    }
+	  else
+	    p_zs(i) = 0.0;
+	}
+    }
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

t2z.h

+/*  t2z.h
+
+    Mark Woolrich & Mark Jenkinson, FMRIB Image Analysis Group
+
+    Copyright (C) 1999-2000 University of Oxford  */
+
+/*  CCOPYRIGHT  */
+
+#if !defined(__t2z_h)
+#define __t2z_h
+  
+#include <iostream>
+#include <fstream>
+#include "newmatap.h"
+#include "newmatio.h"
+#include "base2z.h"
+
+using namespace NEWMAT;
+
+namespace MISCMATHS {
+
+  class T2z : public Base2z
+    {
+    public:
+      static T2z& getInstance();
+      ~T2z() { delete t2z; }
+      
+      float convert(float t, int dof);
+      float converttologp(float t, int dof);
+   
+      static void ComputePs(const ColumnVector& p_vars, const ColumnVector& p_cbs, int p_dof, ColumnVector& p_ps);
+      static void ComputeZStats(const ColumnVector& p_vars, const ColumnVector& p_cbs, int p_dof, ColumnVector& p_zs);
+      static void ComputeZStats(const ColumnVector& p_vars, const ColumnVector& p_cbs, const ColumnVector& p_dof, ColumnVector& p_zs);
+
+    private:
+      T2z() : Base2z()
+	{}
+      
+      const T2z& operator=(T2z&);
+      T2z(T2z&);
+      
+      bool issmalllogp(float logp);
+      bool islarget(float t, int dof, float &logp);
+      float larget2logp(float t, int dof);
+
+      static T2z* t2z;
+
+    };
+ 
+  inline T2z& T2z::getInstance(){
+    if(t2z == NULL)
+      t2z = new T2z();
+  
+    return *t2z;
+  }
+
+}
+
+#endif
+