1
2
3
4
5
6
7
8
9
10
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
fn split_volume_sample(
a: f32,
c: f32, // opacity and colour of original sample
zf: f32,
zb: f32, // front and back of original sample
z: f32, // position of split
) -> ((f32, f32), (f32, f32)) {
// Given a volume sample whose front and back are at depths zf and zb
// respectively, split the sample at depth z. Return the opacities
// and colors of the two parts that result from the split.
//
// The code below is written to avoid excessive rounding errors when
// the opacity of the original sample is very small:
//
// The straightforward computation of the opacity of either part
// requires evaluating an expression of the form
//
// 1.0 - (1.0 - a).pow(x)
//
// However, if a is very small, then 1-a evaluates to 1.0 exactly,
// and the entire expression evaluates to 0.0.
//
// We can avoid this by rewriting the expression as
//
// 1.0 - (x * (1.0 - a).log()).exp()
//
// and replacing the call to log() with a call to the method ln_1p(),
// which computes the logarithm of 1+x without attempting to evaluate
// the expression 1+x when x is very small.
//
// Now we have
//
// 1.0 - (x * (-a).ln_1p()).exp()
//
// However, if a is very small then the call to exp() returns 1.0, and
// the overall expression still evaluates to 0.0. We can avoid that
// by replacing the call to exp() with a call to expm1():
//
// -(x * (-a).ln_1p()).exp_m1()
//
// x.exp_m1() computes exp(x) - 1 in such a way that the result is accurate
// even if x is very small.
//
assert!(zb > zf && z >= zf && z <= zb);
let a = a.clamp(0.0, 1.0);
if a == 1.0f32 {
((1.0f32, c), (1.0f32, c))
} else {
let xf = (z - zf) / (zb - zf);
let xb = (zb - z) / (zb - zf);
if a > f32::MIN_POSITIVE {
//let af = -expm1 (xf * log1p (-a));
let af = -((-a).ln_1p() * xf).exp_m1();
let cf = (af / a) * c;
//ab = -expm1 (xb * log1p (-a));
let ab = -((-a).ln_1p() * xb).exp_m1();
let cb = (ab / a) * c;
((af, cf), (ab, cb))
} else {
((a * xf, c * xf), (a * xb, c * xb))
}
}
}
fn main() {
let a = 0.5f32;
let c = 1.0f32;
let zf = 0.0f32;
let zb = 1.0f32;
assert_eq!(
split_volume_sample(a, c, zf, zb, 0.5),
((0.29289323, 0.58578646), (0.29289323, 0.58578646))
);
assert_eq!(
split_volume_sample(a, c, zf, zb, 1.0e-7),
(
(0.000000069314716, 0.00000013862943),
(0.49999997, 0.99999994)
)
);
}