warped_sampling/
lib.rs

1#[derive(Debug, Copy, Clone)]
2struct Box2f {
3    min: (f64, f64),
4    max: (f64, f64),
5}
6
7fn lerp(a: f64, b: f64, t: f64) -> f64 {
8    a + t * (b - a)
9}
10
11fn lerp_factor(a: f64, b: f64, m: f64) -> f64 {
12    (m - a) / (b - a)
13}
14
15fn set_box(i: usize, edge: [[f64; 2]; 3]) -> Box2f {
16    let ox = i & 1;
17    let oy = (i & 2) >> 1;
18    Box2f {
19        min: (edge[0 + ox][0], edge[0 + oy][1]),
20        max: (edge[1 + ox][0], edge[1 + oy][1]),
21    }
22}
23
24fn warp_a_point(point: (f64, f64), texture_region: &Box2f, warp: &Box2f) -> (f64, f64) {
25    let xf = lerp_factor(warp.min.0, warp.max.0, point.0);
26    let yf = lerp_factor(warp.min.1, warp.max.1, point.1);
27
28    (
29        lerp(texture_region.min.0, texture_region.max.0, xf),
30        lerp(texture_region.min.1, texture_region.max.1, yf),
31    )
32}
33
34pub fn warp(mipmaps: &Vec<&Vec<Vec<f64>>>, points: &[(f64, f64)], warped_points: &mut Vec<(f64, f64)>) {
35    warp_recurse(
36        &mipmaps,
37        0,
38        0,
39        0,
40        2,
41        1.0,
42        &Box2f {
43            min: (0.0, 0.0),
44            max: (1.0, 1.0),
45        },
46        &points,
47        warped_points,
48    );
49}
50
51fn warp_recurse(
52    mipmaps: &Vec<&Vec<Vec<f64>>>,
53    level: usize,
54    x: usize,
55    y: usize,
56    w: usize,
57    scale: f64,
58    b: &Box2f,
59    points: &[(f64, f64)],
60    warped_points: &mut Vec<(f64, f64)>,
61) {
62    if level >= mipmaps.len() {
63        let inv_size = 1.0 / ((1 << level + 1) as f64);
64
65        let texture_box = Box2f {
66            min: (x as f64 * inv_size, y as f64 * inv_size),
67            max: ((x as f64 + 2.0) * inv_size, (y as f64 + 2.0) * inv_size),
68        };
69
70        for (x, y) in points {
71            warped_points.push(warp_a_point((*x, *y), &texture_box, b))
72        }
73
74        return;
75    }
76
77    let cur_mip = &mipmaps[level];
78    let pdfs = [
79        scale * cur_mip[y][x],
80        scale * cur_mip[y][x + 1],
81        scale * cur_mip[y + 1][x],
82        scale * cur_mip[y + 1][x + 1],
83    ];
84
85    let sum = 1.0 / (pdfs[0] + pdfs[1] + pdfs[2] + pdfs[3]);
86    let mid_1 = sum * (pdfs[0] + pdfs[1]);
87    let mid_0 = sum * pdfs[0] / mid_1;
88
89    let mut edge = [
90        [b.min.0, b.min.1],
91        [lerp(b.min.0, b.max.0, mid_0), lerp(b.min.1, b.max.1, mid_1)],
92        [b.max.0, b.max.1],
93    ];
94
95    let b0 = set_box(0, edge);
96    let b1 = set_box(1, edge);
97
98    let mid_2 = sum * pdfs[2] / (1.0 - mid_1);
99    edge[1][0] = lerp(b.min.0, b.max.0, mid_2);
100
101    let b2 = set_box(2, edge);
102    let b3 = set_box(3, edge);
103
104    let boxes = [b0, b1, b2, b3];
105
106    let mut point_cache = [vec![], vec![], vec![], vec![]];
107
108    for (x, y) in points {
109        let mut offset = 0;
110        offset += if *y > boxes[offset].max.1 { 2 } else { 0 };
111        offset += if *x > boxes[offset].max.0 { 1 } else { 0 };
112        point_cache[offset].push((*x, *y));
113    }
114
115    let nl = level + 1;
116    let nx = x << 1;
117    let ny = y << 1;
118    let nw = w << 1;
119
120    for idx in 0..4 {
121        let ox = (idx & 1) << 1;
122        let oy = idx & 2;
123        if point_cache[idx].len() > 0 {
124            warp_recurse(
125                mipmaps,
126                nl,
127                nx + ox,
128                ny + oy,
129                nw,
130                pdfs[idx],
131                &boxes[idx],
132                &point_cache[idx],
133                warped_points,
134            )
135        }
136    }
137}