1use std::f64::consts::PI;
2
3pub const SAMPLES: usize = 81;
4const LOW_NM: f64 = 380.0;
5const STEP_NM: f64 = 5.0;
6
7pub const WEIGHTS: [[f32; 3]; 81] = [
8 [5.039987e-05, 0.0, 0.00035778867],
9 [6.775547e-05, 0.0, 0.0005768482],
10 [0.00013374595, 0.0, 0.001105031],
11 [0.000244155, -8.0429134e-05, 0.001990094],
12 [0.0004620315, -0.00031065758, 0.0037327078],
13 [0.00075178174, -0.00079914776, 0.0060585286],
14 [0.0013901147, -0.0016784245, 0.01140122],
15 [0.0024506692, -0.0030452313, 0.020409657],
16 [0.0041841697, -0.00489408, 0.035484146],
17 [0.006491386, -0.0070713423, 0.05709876],
18 [0.008225769, -0.009268977, 0.07611217],
19 [0.00893323, -0.0110506825, 0.08910742],
20 [0.008651618, -0.011844616, 0.09585401],
21 [0.0075308583, -0.0111710355, 0.09771433],
22 [0.00574333, -0.009713182, 0.0970225],
23 [0.0034815131, -0.0075853216, 0.09534322],
24 [0.00069844024, -0.0048573813, 0.091060266],
25 [-0.0036886954, -0.0017018986, 0.083120264],
26 [-0.0063273897, 0.0017186198, 0.06969027],
27 [-0.008724144, 0.005763426, 0.05593676],
28 [-0.011347991, 0.009961927, 0.04306675],
29 [-0.014388451, 0.014161423, 0.031928584],
30 [-0.017747117, 0.018677758, 0.023228362],
31 [-0.021200482, 0.023959707, 0.016556388],
32 [-0.024560854, 0.030223139, 0.011426527],
33 [-0.027776832, 0.03802389, 0.00728438],
34 [-0.030845057, 0.046441782, 0.0033540253],
35 [-0.03346058, 0.055137716, 0.0],
36 [-0.03508027, 0.06285677, -0.0025753065],
37 [-0.035013195, 0.068301246, -0.004734569],
38 [-0.032630157, 0.071965545, -0.006312869],
39 [-0.027983842, 0.07395785, -0.00737929],
40 [-0.02177047, 0.07446333, -0.008106243],
41 [-0.0140442625, 0.07356939, -0.008540686],
42 [-0.0049150633, 0.07139566, -0.008724192],
43 [0.0046487506, 0.06808836, -0.008696091],
44 [0.0153741175, 0.063680954, -0.008493348],
45 [0.026958833, 0.058147695, -0.008150433],
46 [0.03911687, 0.051680956, -0.007698217],
47 [0.05145111, 0.04444711, -0.0071432022],
48 [0.063475594, 0.036714625, -0.006503886],
49 [0.07455033, 0.028763412, -0.005807879],
50 [0.08411886, 0.020992132, -0.0050818375],
51 [0.091669224, 0.0137893725, -0.00435086],
52 [0.09618923, 0.007610871, -0.0036377995],
53 [0.097940184, 0.002461353, -0.0029763517],
54 [0.09633994, -0.0015017615, -0.0023928043],
55 [0.0919427, -0.004451467, -0.0018944021],
56 [0.08495084, -0.0062283436, -0.0014817494],
57 [0.07555293, -0.0070042815, -0.0011499332],
58 [0.06515909, -0.006994329, -0.00089005305],
59 [0.05535661, -0.0064279693, -0.00069090264],
60 [0.046014983, -0.005523557, -0.0005405573],
61 [0.037232462, -0.004469, -0.00042767118],
62 [0.029350668, -0.0034102271, -0.0003423706],
63 [0.022698376, -0.0024474342, -0.00027671963],
64 [0.017146055, -0.001637843, -0.00022480598],
65 [0.012616245, -0.0010029421, -0.00018254654],
66 [0.00910757, -0.0005379737, -0.00014732922],
67 [0.0066325837, -0.00022167167, -0.00011759913],
68 [0.004884865, -2.4797864e-05, -9.247229e-05],
69 [0.0034370946, 0.0, -7.142631e-05],
70 [0.0023721277, 0.0, -5.40864e-05],
71 [0.0016515272, 0.0, -4.010171e-05],
72 [0.0011923393, 0.0, -2.9092696e-05],
73 [0.00084797293, 0.0, -2.064557e-05],
74 [0.000605779, 0.0, -1.4331459e-05],
75 [0.0004273203, 0.0, -9.733269e-06],
76 [0.00030592916, 0.0, -6.469539e-06],
77 [0.00021032631, 0.0, -4.210301e-06],
78 [0.00014668911, 0.0, -2.683982e-06],
79 [0.0001021725, 0.0, -1.6768181e-06],
80 [7.628618e-05, 0.0, -1.027186e-06],
81 [5.108625e-05, 0.0, -6.1728224e-07],
82 [3.186762e-05, 0.0, -3.6408065e-07],
83 [1.9218627e-05, 0.0, -2.1085688e-07],
84 [1.9218627e-05, 0.0, -1.1996055e-07],
85 [1.2648995e-05, 0.0, 2.5954927e-07],
86 [1.2648995e-05, 0.0, 2.5954927e-07],
87 [1.2648995e-05, 0.0, 2.5954927e-07],
88 [-5.4965028e-08, 0.0, -1.0592136e-08],
89];
90
91fn fresnel(cos_i: f64, from: f64, to: f64) -> Option<(f64, f64, f64)> {
92 let sin2 = (from / to) * (from / to) * (1.0 - cos_i * cos_i);
93 if sin2 >= 1.0 {
94 return None;
95 }
96 let cos_t = (1.0 - sin2).sqrt();
97 let s = (from * cos_i - to * cos_t) / (from * cos_i + to * cos_t);
98 let p = (to * cos_i - from * cos_t) / (to * cos_i + from * cos_t);
99 Some((s, p, cos_t))
100}
101
102pub fn film_reflectance(
103 cosine: f32,
104 thickness_nm: f32,
105 outer_ior: f32,
106 film_ior: f32,
107 inner_ior: f32,
108) -> [f32; 3] {
109 let (outer, film, inner) = (outer_ior as f64, film_ior as f64, inner_ior as f64);
110 let Some(first) = fresnel((cosine as f64).clamp(0.0, 1.0), outer, film) else {
111 return [1.0; 3];
112 };
113 let Some(second) = fresnel(first.2, film, inner) else {
114 return [1.0; 3];
115 };
116 let path = 4.0 * PI * film * thickness_nm as f64 * first.2;
117 let r12 = [first.0, first.1];
118 let r23 = [second.0, second.1];
119 let a = [
120 r12[0] * r12[0] + r23[0] * r23[0],
121 r12[1] * r12[1] + r23[1] * r23[1],
122 ];
123 let b = [2.0 * r12[0] * r23[0], 2.0 * r12[1] * r23[1]];
124 let e = [
125 1.0 + r12[0] * r12[0] * r23[0] * r23[0],
126 1.0 + r12[1] * r12[1] * r23[1] * r23[1],
127 ];
128 let mut total = [0.0; 3];
129 for (k, weight) in WEIGHTS.iter().enumerate() {
130 let c = (path / (LOW_NM + STEP_NM * k as f64)).cos();
131 let reflect =
132 0.5 * ((a[0] + b[0] * c) / (e[0] + b[0] * c) + (a[1] + b[1] * c) / (e[1] + b[1] * c));
133 for channel in 0..3 {
134 total[channel] += weight[channel] as f64 * reflect;
135 }
136 }
137 total.map(|v| v.clamp(0.0, 1.0) as f32)
138}
139
140fn product(a: (f64, f64), b: (f64, f64)) -> (f64, f64) {
141 (a.0 * b.0 - a.1 * b.1, a.0 * b.1 + a.1 * b.0)
142}
143
144fn quotient(a: (f64, f64), b: (f64, f64)) -> (f64, f64) {
145 let size = b.0 * b.0 + b.1 * b.1;
146 (
147 (a.0 * b.0 + a.1 * b.1) / size,
148 (a.1 * b.0 - a.0 * b.1) / size,
149 )
150}
151
152fn root(a: (f64, f64)) -> (f64, f64) {
153 let size = a.0.hypot(a.1);
154 let im = (0.5 * (size - a.0)).max(0.0).sqrt();
155 (
156 (0.5 * (size + a.0)).max(0.0).sqrt(),
157 if a.1 < 0.0 { -im } else { im },
158 )
159}
160
161pub fn conductor_reflectance(
162 cosine: f32,
163 thickness_nm: f32,
164 film_ior: f32,
165 f0: [f32; 3],
166) -> [f32; 3] {
167 let film = film_ior.max(1.01) as f64;
168 let Some((rs12, rp12, cos_f)) = fresnel((cosine as f64).clamp(0.0, 1.0), 1.0, film) else {
169 return [1.0; 3];
170 };
171 let path = 4.0 * PI * film * thickness_nm as f64 * cos_f;
172 std::array::from_fn(|channel| {
173 let [n, k] = pfx_materials::conductor_ior(f0[channel]).map(f64::from);
174 let square = product((n, k), (n, k));
175 let u = root((square.0 - film * film * (1.0 - cos_f * cos_f), square.1));
176 let across = film * cos_f;
177 let rs23 = quotient((across - u.0, -u.1), (across + u.0, u.1));
178 let tilted = (square.0 * cos_f, square.1 * cos_f);
179 let rp23 = quotient(
180 (tilted.0 - film * u.0, tilted.1 - film * u.1),
181 (tilted.0 + film * u.0, tilted.1 + film * u.1),
182 );
183 let mut total = 0.0;
184 for (step, weight) in WEIGHTS.iter().enumerate() {
185 let phase = path / (LOW_NM + STEP_NM * step as f64);
186 let z = (phase.cos(), phase.sin());
187 let mut reflect = 0.0;
188 for (r12, r23) in [(rs12, rs23), (rp12, rp23)] {
189 let q = product(r23, z);
190 let top = (r12 + q.0).powi(2) + q.1.powi(2);
191 let bottom = (1.0 + r12 * q.0).powi(2) + (r12 * q.1).powi(2);
192 reflect += 0.5 * top / bottom;
193 }
194 total += weight[channel] as f64 * reflect;
195 }
196 total.clamp(0.0, 1.0) as f32
197 })
198}
199
200pub fn surface_reflectance(
201 cosine: f32,
202 thickness_nm: f32,
203 film_ior: f32,
204 ior: f32,
205 entering: bool,
206) -> [f32; 3] {
207 let film = film_ior.max(1.01);
208 if entering {
209 film_reflectance(cosine, thickness_nm, 1.0, film, ior)
210 } else {
211 film_reflectance(cosine, thickness_nm, 1.0, film / ior, 1.0 / ior)
212 }
213}
214
215pub fn srgb_byte(linear: f32) -> f32 {
216 let v = linear.clamp(0.0, 1.0);
217 255.0 * pfx_materials::encode_channel(v)
218}
219
220#[cfg(test)]
221mod tests {
222 use super::*;
223
224 #[test]
225 fn a_flat_spectrum_stays_white() {
226 for channel in 0..3 {
227 let sum: f64 = WEIGHTS.iter().map(|row| row[channel] as f64).sum();
228 assert!((sum - 1.0).abs() < 1e-6);
229 }
230 }
231
232 #[test]
233 fn the_shader_table_is_the_same_table() {
234 let source = include_str!("trace.wgsl");
235 let start = source.find("const CF_WEIGHTS").unwrap();
236 let end = start + source[start..].find(");").unwrap();
237 let numbers: Vec<f64> = source[start..end]
238 .split("vec3f(")
239 .skip(1)
240 .flat_map(|row| {
241 row.split(')')
242 .next()
243 .unwrap()
244 .split(',')
245 .map(|n| n.trim().parse::<f64>().unwrap())
246 .collect::<Vec<_>>()
247 })
248 .collect();
249 assert_eq!(numbers.len(), SAMPLES * 3);
250 for (k, row) in WEIGHTS.iter().enumerate() {
251 for channel in 0..3 {
252 assert_eq!(numbers[k * 3 + channel] as f32, row[channel]);
253 }
254 }
255 }
256
257 #[test]
258 fn no_film_is_the_plain_fresnel() {
259 let thin = film_reflectance(1.0, 0.001, 1.0, 1.9, 1.4);
260 for value in thin {
261 assert!((value - 0.0278).abs() < 1e-3, "{value}");
262 }
263 }
264
265 #[test]
266 fn leaving_reads_the_film_thinner_by_the_ior_as_cycles_does() {
267 for cosine in [1.0, 0.8, 0.5, 0.2] {
268 for thickness in [18.0, 90.0, 180.0, 400.0] {
269 let ours = surface_reflectance(cosine, thickness, 1.9, 1.4, false);
270 let thinner = film_reflectance(cosine, thickness / 1.4, 1.4, 1.9, 1.0);
271 for channel in 0..3 {
272 assert!(
273 (ours[channel] - thinner[channel]).abs() < 1e-5,
274 "{ours:?} {thinner:?}"
275 );
276 }
277 }
278 }
279 let physical = film_reflectance(1.0, 180.0, 1.4, 1.9, 1.0);
280 let leaving = surface_reflectance(1.0, 180.0, 1.9, 1.4, false);
281 assert!((0..3).any(|c| (physical[c] - leaving[c]).abs() > 0.01));
282 assert_eq!(
283 surface_reflectance(0.7, 180.0, 1.9, 1.4, true),
284 film_reflectance(0.7, 180.0, 1.0, 1.9, 1.4)
285 );
286 }
287
288 #[test]
289 fn a_bare_conductor_is_its_complex_fresnel() {
290 let f0 = [0.95, 0.64, 0.54];
291 let bare = conductor_reflectance(1.0, 0.0, 1.6, f0);
292 for channel in 0..3 {
293 assert!((bare[channel] - f0[channel]).abs() < 1e-5, "{bare:?}");
294 }
295 let thin = conductor_reflectance(1.0, 120.0, 1.6, f0);
296 assert!((0..3).any(|c| (thin[c] - f0[c]).abs() > 0.02), "{thin:?}");
297 }
298}