Newer
Older
/* base2z.cc
Mark Woolrich & Mark Jenkinson, FMRIB Image Analysis Group
Copyright (C) 1999-2000 University of Oxford */
/* CCOPYRIGHT */
#include "base2z.h"
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
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;
}
}