Newer
Older
/* base2z.cc
Mark Woolrich & Mark Jenkinson, FMRIB Image Analysis Group
Copyright (C) 1999-2000 University of Oxford */
/* CCOPYRIGHT */
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
101
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;
}
}