1// Copyright (c) 2022-2022, The rav1e contributors. All rights reserved
2//
3// This source code is subject to the terms of the BSD 2 Clause License and
4// the Alliance for Open Media Patent License 1.0. If the BSD 2 Clause License
5// was not distributed with this source code in the LICENSE file, you can
6// obtain it at www.aomedia.org/license/software. If the Alliance for Open
7// Media Patent License 1.0 was not distributed with this source code in the
8// PATENTS file, you can obtain it at www.aomedia.org/license/patent.
9
10// The original work for this formula was implmented in aomenc, and this is
11// an adaptation of that work:
12// https://aomedia.googlesource.com/aom/+/refs/heads/main/examples/photon_noise_table.c
13
14// This implementation creates a film grain table, for use in stills and videos,
15// representing the noise that one would get by shooting with a digital camera
16// at a given light level. Much of the noise in digital images is photon shot
17// noise, which is due to the characteristics of photon arrival and grows in
18// standard deviation as the square root of the expected number of photons
19// captured.
20// https://www.photonstophotos.net/Emil%20Martinec/noise.html#shotnoise
21//
22// The proxy used by this implementation for the amount of light captured
23// is the ISO value such that the focal plane exposure at the time of capture
24// would have been mapped by a 35mm camera to the output lightness observed
25// in the image. That is, if one were to shoot on a 35mm camera (36×24mm sensor)
26// at the nominal exposure for that ISO setting, the resulting image should
27// contain noise of the same order of magnitude as generated by this
28// implementation.
29//
30// The (mostly) square-root relationship between light intensity and noise
31// amplitude holds in linear light, but AV1 streams are most often encoded
32// non-linearly, and the film grain is applied to those non-linear values.
33// Therefore, this implementation must account for the non-linearity, and this
34// is controlled by the transfer function parameter, which specifies the tone
35// response curve that will be used when encoding the actual image. The default
36// for this implementation is BT.1886, which is approximately similar to an
37// encoding gamma of 1/2.8 (i.e. a decoding gamma of 2.8) though not quite
38// identical.
39//
40// As alluded to above, the implementation assumes that the image is taken from
41// the entirety of a 36×24mm (“35mm format”) sensor. If that assumption does not
42// hold, then a “35mm-equivalent ISO value” that can be passed to the
43// implementation can be obtained by multiplying the true ISO value by the ratio
44// of 36×24mm to the area that was actually used. For formats that approximately
45// share the same aspect ratio, this is often expressed as the square of the
46// “equivalence ratio” which is the ratio of their diagonals. For example, APS-C
47// (often ~24×16mm) is said to have an equivalence ratio of 1.5 relative to the
48// 35mm format, and therefore ISO 1000 on APS-C and ISO 1000×1.5² = 2250 on 35mm
49// produce an image of the same lightness from the same amount of light spread
50// onto their respective surface areas (resulting in different focal plane
51// exposures), and those images will thus have similar amounts of noise if the
52// cameras are of similar technology. https://doi.org/10.1117/1.OE.57.11.110801
53//
54// The implementation needs to know the resolution of the images to which its
55// grain tables will be applied so that it can know how the light on the sensor
56// was shared between its pixels. As a general rule, while a higher pixel count
57// will lead to more noise per pixel, when the final image is viewed at the same
58// physical size, that noise will tend to “average out” to the same amount over
59// a given area, since there will be more pixels in it which, in aggregate, will
60// have received essentially as much light. Put differently, the amount of noise
61// depends on the scale at which it is measured, and the decision for this
62// implementation was to make that scale relative to the image instead of its
63// constituent samples. For more on this, see:
64//
65// https://www.photonstophotos.net/Emil%20Martinec/noise-p3.html#pixelsize
66// https://www.dpreview.com/articles/5365920428/the-effect-of-pixel-and-sensor-sizes-on-noise/2
67// https://www.dpreview.com/videos/7940373140/dpreview-tv-why-lower-resolution-sensors-are-not-better-in-low-light
68
69use std::{
70 fs::File,
71 io::{BufWriter, Write},
72 path::Path,
73};
74
75use arrayvec::ArrayVec;
76
77use crate::{GrainTableSegment, ScalingPoints, DEFAULT_GRAIN_SEED, NUM_Y_POINTS};
78
79const PQ_M1: f32 = 2610. / 16384.;
80const PQ_M2: f32 = 128. * 2523. / 4096.;
81const PQ_C1: f32 = 3424. / 4096.;
82const PQ_C2: f32 = 32. * 2413. / 4096.;
83const PQ_C3: f32 = 32. * 2392. / 4096.;
84
85const BT1886_WHITEPOINT: f32 = 203.;
86const BT1886_BLACKPOINT: f32 = 0.1;
87const BT1886_GAMMA: f32 = 2.4;
88
89// BT.1886 formula from https://en.wikipedia.org/wiki/ITU-R_BT.1886.
90//
91// TODO: the inverses, alpha, and beta should all be constants
92// once floats in const fns are stabilized and `powf` is const.
93// Until then, `inline(always)` gets us close enough.
94
95#[inline(always)]
96fn bt1886_inv_whitepoint() -> f32 {
97 BT1886_WHITEPOINT.powf(1.0 / BT1886_GAMMA)
98}
99
100#[inline(always)]
101fn bt1886_inv_blackpoint() -> f32 {
102 BT1886_BLACKPOINT.powf(1.0 / BT1886_GAMMA)
103}
104
105/// The variable for user gain:
106/// `α = (Lw^(1/λ) - Lb^(1/λ)) ^ λ`
107#[inline(always)]
108fn bt1886_alpha() -> f32 {
109 (bt1886_inv_whitepoint() - bt1886_inv_blackpoint()).powf(BT1886_GAMMA)
110}
111
112/// The variable for user black level lift:
113/// `β = Lb^(1/λ) / (Lw^(1/λ) - Lb^(1/λ))`
114#[inline(always)]
115fn bt1886_beta() -> f32 {
116 bt1886_inv_blackpoint() / (bt1886_inv_whitepoint() - bt1886_inv_blackpoint())
117}
118
119/// Settings and video data defining how to generate the film grain params.
120#[derive(Debug, Clone, Copy)]
121pub struct NoiseGenArgs {
122 pub iso_setting: u32,
123 pub width: u32,
124 pub height: u32,
125 pub transfer_function: TransferFunction,
126 pub chroma_grain: bool,
127 pub random_seed: Option<u16>,
128}
129
130/// Generates a set of photon noise parameters for a segment of video
131/// given a set of `args`.
132#[must_use]
133pub fn generate_photon_noise_params(
134 start_time: u64,
135 end_time: u64,
136 args: NoiseGenArgs,
137) -> GrainTableSegment {
138 GrainTableSegment {
139 start_time,
140 end_time,
141 scaling_points_y: generate_luma_noise_points(args),
142 scaling_points_cb: ArrayVec::new(),
143 scaling_points_cr: ArrayVec::new(),
144 scaling_shift: 8,
145 ar_coeff_lag: 0,
146 ar_coeffs_y: ArrayVec::new(),
147 ar_coeffs_cb: ArrayVec::try_from([0].as_slice())
148 .expect("Cannot fail creation from const array"),
149 ar_coeffs_cr: ArrayVec::try_from([0].as_slice())
150 .expect("Cannot fail creation from const array"),
151 ar_coeff_shift: 6,
152 cb_mult: 0,
153 cb_luma_mult: 0,
154 cb_offset: 0,
155 cr_mult: 0,
156 cr_luma_mult: 0,
157 cr_offset: 0,
158 overlap_flag: true,
159 chroma_scaling_from_luma: args.chroma_grain,
160 grain_scale_shift: 0,
161 random_seed: args.random_seed.unwrap_or(DEFAULT_GRAIN_SEED),
162 }
163}
164
165/// Generates a set of film grain parameters for a segment of video
166/// given a set of `args`.
167///
168/// # Panics
169/// - This is not yet implemented, so it will always panic
170#[must_use]
171#[cfg(feature = "unstable")]
172pub fn generate_film_grain_params(
173 start_time: u64,
174 end_time: u64,
175 args: NoiseGenArgs,
176) -> GrainTableSegment {
177 todo!("SCIENCE");
178 // GrainTableSegment {
179 // start_time,
180 // end_time,
181 // scaling_points_y: generate_luma_noise_points(args),
182 // scaling_points_cb: ArrayVec::new(),
183 // scaling_points_cr: ArrayVec::new(),
184 // scaling_shift: 8,
185 // ar_coeff_lag: 0,
186 // ar_coeffs_y: ArrayVec::new(),
187 // ar_coeffs_cb: ArrayVec::try_from([0].as_slice())
188 // .expect("Cannot fail creation from const array"),
189 // ar_coeffs_cr: ArrayVec::try_from([0].as_slice())
190 // .expect("Cannot fail creation from const array"),
191 // ar_coeff_shift: 6,
192 // cb_mult: 0,
193 // cb_luma_mult: 0,
194 // cb_offset: 0,
195 // cr_mult: 0,
196 // cr_luma_mult: 0,
197 // cr_offset: 0,
198 // overlap_flag: true,
199 // chroma_scaling_from_luma: args.chroma_grain,
200 // grain_scale_shift: 0,
201 // random_seed: args.random_seed.unwrap_or(DEFAULT_GRAIN_SEED),
202 // }
203}
204
205/// Write a set of generated film grain params to a table file,
206/// using the standard film grain table format supported by
207/// aomenc, rav1e, and svt-av1.
208///
209/// # Errors
210///
211/// - If the output file cannot be written to
212pub fn write_grain_table<P: AsRef<Path>>(
213 filename: P,
214 params: &[GrainTableSegment],
215) -> anyhow::Result<()> {
216 let mut file: BufWriter = BufWriter::new(inner:File::create(path:filename)?);
217 writeln!(&mut file, "filmgrn1")?;
218 for segment: &GrainTableSegment in params {
219 write_film_grain_segment(params:segment, &mut file)?;
220 }
221 file.flush()?;
222
223 Ok(())
224}
225
226fn write_film_grain_segment(
227 params: &GrainTableSegment,
228 output: &mut BufWriter<File>,
229) -> anyhow::Result<()> {
230 writeln!(
231 output,
232 "E {} {} 1 {} 1",
233 params.start_time, params.end_time, params.random_seed,
234 )?;
235 writeln!(
236 output,
237 "\tp {} {} {} {} {} {} {} {} {} {} {} {}",
238 params.ar_coeff_lag,
239 params.ar_coeff_shift,
240 params.grain_scale_shift,
241 params.scaling_shift,
242 u8::from(params.chroma_scaling_from_luma),
243 u8::from(params.overlap_flag),
244 params.cb_mult,
245 params.cb_luma_mult,
246 params.cb_offset,
247 params.cr_mult,
248 params.cr_luma_mult,
249 params.cr_offset
250 )?;
251
252 write!(output, "\tsY {} ", params.scaling_points_y.len())?;
253 for point in &params.scaling_points_y {
254 write!(output, " {} {}", point[0], point[1])?;
255 }
256 writeln!(output)?;
257
258 write!(output, "\tsCb {}", params.scaling_points_cb.len())?;
259 for point in &params.scaling_points_cb {
260 write!(output, " {} {}", point[0], point[1])?;
261 }
262 writeln!(output)?;
263
264 write!(output, "\tsCr {}", params.scaling_points_cr.len())?;
265 for point in &params.scaling_points_cr {
266 write!(output, " {} {}", point[0], point[1])?;
267 }
268 writeln!(output)?;
269
270 write!(output, "\tcY")?;
271 for coeff in &params.ar_coeffs_y {
272 write!(output, " {}", *coeff)?;
273 }
274 writeln!(output)?;
275
276 write!(output, "\tcCb")?;
277 for coeff in &params.ar_coeffs_cb {
278 write!(output, " {}", *coeff)?;
279 }
280 writeln!(output)?;
281
282 write!(output, "\tcCr")?;
283 for coeff in &params.ar_coeffs_cr {
284 write!(output, " {}", *coeff)?;
285 }
286 writeln!(output)?;
287
288 Ok(())
289}
290
291#[allow(clippy::upper_case_acronyms)]
292#[derive(Debug, Clone, Copy, PartialEq, Eq)]
293pub enum TransferFunction {
294 /// For SDR content
295 BT1886,
296 /// For HDR content
297 SMPTE2084,
298}
299
300impl TransferFunction {
301 #[must_use]
302 pub fn to_linear(self, x: f32) -> f32 {
303 match self {
304 TransferFunction::BT1886 => {
305 // The screen luminance in cd/m^2:
306 // L = α * (x + β)^λ
307 let luma = bt1886_alpha() * (x + bt1886_beta()).powf(BT1886_GAMMA);
308
309 // Normalize to between 0.0 and 1.0
310 luma / BT1886_WHITEPOINT
311 }
312 TransferFunction::SMPTE2084 => {
313 let pq_pow_inv_m2 = x.powf(1. / PQ_M2);
314 (0_f32.max(pq_pow_inv_m2 - PQ_C1) / PQ_C3.mul_add(-pq_pow_inv_m2, PQ_C2))
315 .powf(1. / PQ_M1)
316 }
317 }
318 }
319
320 #[allow(clippy::wrong_self_convention)]
321 #[must_use]
322 pub fn from_linear(self, x: f32) -> f32 {
323 match self {
324 TransferFunction::BT1886 => {
325 // Scale to a raw cd/m^2 value
326 let luma = x * BT1886_WHITEPOINT;
327
328 // The inverse of the `to_linear` formula:
329 // `(L / α)^(1 / λ) - β = x`
330 (luma / bt1886_alpha()).powf(1.0 / BT1886_GAMMA) - bt1886_beta()
331 }
332 TransferFunction::SMPTE2084 => {
333 if x < f32::EPSILON {
334 return 0.0;
335 }
336 let linear_pow_m1 = x.powf(PQ_M1);
337 (PQ_C2.mul_add(linear_pow_m1, PQ_C1) / PQ_C3.mul_add(linear_pow_m1, 1.)).powf(PQ_M2)
338 }
339 }
340 }
341
342 #[inline(always)]
343 #[must_use]
344 pub fn mid_tone(self) -> f32 {
345 self.to_linear(0.5)
346 }
347}
348
349fn generate_luma_noise_points(args: NoiseGenArgs) -> ScalingPoints {
350 // Assumes a daylight-like spectrum.
351 // https://www.strollswithmydog.com/effective-quantum-efficiency-of-sensor/#:~:text=11%2C260%20photons/um%5E2/lx-s
352 const PHOTONS_PER_SQ_MICRON_PER_LUX_SECOND: f32 = 11260.;
353
354 // Order of magnitude for cameras in the 2010-2020 decade, taking the CFA into
355 // account.
356 const EFFECTIVE_QUANTUM_EFFICIENCY: f32 = 0.2;
357
358 // Also reasonable values for current cameras. The read noise is typically
359 // higher than this at low ISO settings but it matters less there.
360 const PHOTO_RESPONSE_NON_UNIFORMITY: f32 = 0.005;
361 const INPUT_REFERRED_READ_NOISE: f32 = 1.5;
362
363 // Assumes a 35mm sensor (36mm × 24mm).
364 const SENSOR_AREA: f32 = 36_000. * 24_000.;
365
366 // Focal plane exposure for a mid-tone (typically a 18% reflectance card), in
367 // lx·s.
368 let mid_tone_exposure = 10. / args.iso_setting as f32;
369
370 let pixel_area_microns = SENSOR_AREA / (args.width * args.height) as f32;
371
372 let mid_tone_electrons_per_pixel = EFFECTIVE_QUANTUM_EFFICIENCY
373 * PHOTONS_PER_SQ_MICRON_PER_LUX_SECOND
374 * mid_tone_exposure
375 * pixel_area_microns;
376 let max_electrons_per_pixel = mid_tone_electrons_per_pixel / args.transfer_function.mid_tone();
377
378 let mut scaling_points = ScalingPoints::default();
379 for i in 0..NUM_Y_POINTS {
380 let x = i as f32 / (NUM_Y_POINTS as f32 - 1.);
381 let linear = args.transfer_function.to_linear(x);
382 let electrons_per_pixel = max_electrons_per_pixel * linear;
383
384 // Quadrature sum of the relevant sources of noise, in electrons rms. Photon
385 // shot noise is sqrt(electrons) so we can skip the square root and the
386 // squaring.
387 // https://en.wikipedia.org/wiki/Addition_in_quadrature
388 // https://doi.org/10.1117/3.725073
389 let noise_in_electrons = (PHOTO_RESPONSE_NON_UNIFORMITY
390 * PHOTO_RESPONSE_NON_UNIFORMITY
391 * electrons_per_pixel)
392 .mul_add(
393 electrons_per_pixel,
394 INPUT_REFERRED_READ_NOISE.mul_add(INPUT_REFERRED_READ_NOISE, electrons_per_pixel),
395 )
396 .sqrt();
397 let linear_noise = noise_in_electrons / max_electrons_per_pixel;
398 let linear_range_start = 0_f32.max(2.0f32.mul_add(-linear_noise, linear));
399 let linear_range_end = 1_f32.min(2_f32.mul_add(linear_noise, linear));
400 let tf_slope = (args.transfer_function.from_linear(linear_range_end)
401 - args.transfer_function.from_linear(linear_range_start))
402 / (linear_range_end - linear_range_start);
403 let encoded_noise = linear_noise * tf_slope;
404
405 let x = (255. * x).round() as u8;
406 let encoded_noise = 255_f32.min((255. * 7.88 * encoded_noise).round()) as u8;
407
408 scaling_points.push([x, encoded_noise]);
409 }
410
411 scaling_points
412}
413
414#[cfg(test)]
415mod tests {
416 use quickcheck::TestResult;
417 use quickcheck_macros::quickcheck;
418
419 use super::*;
420
421 #[quickcheck]
422 fn bt1886_to_linear_within_range(x: f32) -> TestResult {
423 if !(0.0..=1.0).contains(&x) || x.is_nan() {
424 return TestResult::discard();
425 }
426
427 let tx = TransferFunction::BT1886;
428 let res = tx.to_linear(x);
429 TestResult::from_bool((0.0..=1.0).contains(&res))
430 }
431
432 #[quickcheck]
433 fn bt1886_to_linear_reverts_correctly(x: f32) -> TestResult {
434 if !(0.0..=1.0).contains(&x) || x.is_nan() {
435 return TestResult::discard();
436 }
437
438 let tx = TransferFunction::BT1886;
439 let res = tx.to_linear(x);
440 let res = tx.from_linear(res);
441 TestResult::from_bool((x - res).abs() < f32::EPSILON)
442 }
443
444 #[quickcheck]
445 fn smpte2084_to_linear_within_range(x: f32) -> TestResult {
446 if !(0.0..=1.0).contains(&x) || x.is_nan() {
447 return TestResult::discard();
448 }
449
450 let tx = TransferFunction::SMPTE2084;
451 let res = tx.to_linear(x);
452 TestResult::from_bool((0.0..=1.0).contains(&res))
453 }
454
455 #[quickcheck]
456 fn smpte2084_to_linear_reverts_correctly(x: f32) -> TestResult {
457 if !(0.0..=1.0).contains(&x) || x.is_nan() {
458 return TestResult::discard();
459 }
460
461 let tx = TransferFunction::SMPTE2084;
462 let res = tx.to_linear(x);
463 let res = tx.from_linear(res);
464 TestResult::from_bool((x - res).abs() < f32::EPSILON)
465 }
466}
467