From 1dae8c9fc138d36c5d56d0b7bf61df1663ee7fac Mon Sep 17 00:00:00 2001 From: Nathan Vegdahl Date: Fri, 22 Apr 2022 13:17:56 -0700 Subject: [PATCH] Major performance improvements to transfer function formula estimation. It also now ensures that the end meets exactly where it does in the LUT. --- Cargo.lock | 40 ---------------- Cargo.toml | 1 - src/linear_log.rs | 27 +++++++++++ src/main.rs | 14 ++++-- src/optimize_log.rs | 109 ++++++++++++++++++++++++++++++++------------ 5 files changed, 119 insertions(+), 72 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index 440d4c6..d0cddfa 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -211,7 +211,6 @@ dependencies = [ "clap", "colorbox", "exr", - "simplers_optimization", ] [[package]] @@ -229,15 +228,6 @@ dependencies = [ "getrandom", ] -[[package]] -name = "num-traits" -version = "0.2.14" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "9a64b1ec5cda2586e284722486d802acf1f7dbdc623e2bfc57e65ca1cd099290" -dependencies = [ - "autocfg", -] - [[package]] name = "num_cpus" version = "1.13.1" @@ -248,15 +238,6 @@ dependencies = [ "libc", ] -[[package]] -name = "ordered-float" -version = "2.10.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "7940cf2ca942593318d07fcf2596cdca60a85c9e7fab408a5e21a4f9dcd40d87" -dependencies = [ - "num-traits", -] - [[package]] name = "os_str_bytes" version = "6.0.0" @@ -286,16 +267,6 @@ dependencies = [ "syn", ] -[[package]] -name = "priority-queue" -version = "1.2.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "00ba480ac08d3cfc40dea10fd466fd2c14dee3ea6fc7873bc4079eda2727caf0" -dependencies = [ - "autocfg", - "indexmap", -] - [[package]] name = "proc-macro2" version = "1.0.37" @@ -320,17 +291,6 @@ version = "1.1.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "d29ab0c6d3fc0ee92fe66e2d99f700eab17a8d57d1c1d3b748380fb20baa78cd" -[[package]] -name = "simplers_optimization" -version = "0.4.3" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "2cd97912bb2a16575a2c632c2a2f2bac8a527827ceaddd73e0ccc12d86adec43" -dependencies = [ - "num-traits", - "ordered-float", - "priority-queue", -] - [[package]] name = "smallvec" version = "1.8.0" diff --git a/Cargo.toml b/Cargo.toml index 119a28d..fd38698 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -7,4 +7,3 @@ edition = "2021" exr = "1.4.1" clap = { version = "3.1.8", default-features = false, features=["std"] } colorbox = { git = "https://github.com/cessen/colorbox", branch = "master" } -simplers_optimization = "0.4.3" diff --git a/src/linear_log.rs b/src/linear_log.rs index 300875b..1701dbb 100644 --- a/src/linear_log.rs +++ b/src/linear_log.rs @@ -38,6 +38,33 @@ pub fn log_to_linear(x: f64, line_offset: f64, slope: f64, log_offset: f64, base } } +// Find the `log_offset` needed to put x=end at y=1.0 in the linear_to_log function. +pub fn find_log_offset_for_end(end: f64, line_offset: f64, slope: f64, base: f64) -> f64 { + let mut offset_up = 10.0; + let mut offset_down = -10.0; + + for _ in 0..54 { + let log_offset = (offset_up + offset_down) * 0.5; + if linear_to_log(end, line_offset, slope, log_offset, base) > 1.0 { + offset_up = log_offset; + } else { + offset_down = log_offset; + } + } + + offset_up +} + +// Transition point between log and linear. +// +// Returned as (linear, non-linear). +pub fn transition_point(line_offset: f64, slope: f64, log_offset: f64, base: f64) -> (f64, f64) { + let transition = 1.0 / (slope * base.ln()); + let k = transition + log_offset; + + (k, (k - line_offset) * slope) +} + //------------------------------------------------------------- /// Generates Rust code for a linear-to-log transfer function with the diff --git a/src/main.rs b/src/main.rs index 317b4f9..76a12d0 100644 --- a/src/main.rs +++ b/src/main.rs @@ -62,8 +62,18 @@ fn main() { image }; - // Build the LUT. + // Fetch the transfer function LUT. let gray = &mut input_image[test_image::gray_idx(0)..test_image::gray_idx(GRADIENT_LEN)]; + + // Attempt to find an analytic log-linear function that matches + // the transfer function. + let full_lut: Vec = gray + .iter() + .map(|rgb| ((rgb[0] as f64 + rgb[1] as f64 + rgb[2] as f64) / 3.0) as f32) + .collect(); + optimize_log::find_parameters(&full_lut); + + // Build the LUT for export. let mut prev = gray[0]; for rgb in gray.iter_mut() { // Ensure montonicity. @@ -89,8 +99,6 @@ fn main() { gray_b.push(rgb[2]); } - optimize_log::find_parameters(&gray_r); - // Write the LUT. colorbox::formats::cube::write_1d( BufWriter::new(File::create("test.cube").unwrap()), diff --git a/src/optimize_log.rs b/src/optimize_log.rs index bf6f977..a9e6158 100644 --- a/src/optimize_log.rs +++ b/src/optimize_log.rs @@ -5,41 +5,48 @@ pub fn find_parameters(lut: &[f32]) { // Compute the stuff that we can without estimation. let offset = lut[0] as f64; - let slope = lin_norm / (lut[1] as f64 - lut[0] as f64); + let end = lut[lut.len() - 1] as f64; + let slope = { + // We take the difference of points near zero for increased accuracy. + let (i, _) = lut.iter().enumerate().find(|(_, y)| **y > 0.0).unwrap(); + lin_norm / (lut[i] as f64 - lut[i - 1] as f64) + }; - // Select a range of points from the lookup table to fit to. - let idxs: Vec<_> = (0..lut.len()).step_by(lut.len() / 256).collect(); - let coords: Vec<(f64, f64)> = idxs + // Collect LUT points as (x, y) coordinates. + let coords: Vec<(f64, f64)> = lut .iter() - .map(|i| (*i as f64 * lin_norm, lut[*i] as f64)) + .enumerate() + .step_by(lut.len() / 256) + .map(|(i, y)| (i as f64 * lin_norm, *y as f64)) .collect(); // Do the fitting. - let f = |v: &[f64]| { - let mut avg_sqr_err = 0.0f64; - for (x, y) in coords.iter().copied() { - let e = (log_to_lin(x, offset, slope, v[0], v[1]) - y).abs() / y.abs(); - avg_sqr_err += e * e; - } - let last_y = lut[lut.len() - 1] as f64; - let e = (log_to_lin(1.0, offset, slope, v[0], v[1]) - last_y).abs() / last_y.abs(); - avg_sqr_err += e * e; - avg_sqr_err - }; - let input_interval = vec![(-0.2, 0.2), (1.1, 1000.0)]; - let (_, params) = simplers_optimization::Optimizer::minimize(&f, &input_interval, 1000000); + let base = optimize( + |v: f64| { + let log_offset = crate::linear_log::find_log_offset_for_end(end, offset, slope, v); + let mut sqr_err = 0.0f64; + for (x, y) in coords.iter().copied() { + let e = (log_to_lin(x, offset, slope, log_offset, v) - y).abs() / y.abs(); + sqr_err += e * e; + } + sqr_err + }, + [1.5, 10000000.0], + 100, + ); + let log_offset = crate::linear_log::find_log_offset_for_end(end, offset, slope, base); + let transition = crate::linear_log::transition_point(offset, slope, log_offset, base); // Calculate the error of our model. let mut max_err = 0.0f64; let mut avg_err = 0.0f64; let mut avg_samples = 0usize; - for (i, y) in lut.iter().map(|y| *y as f64).enumerate() { - // We only record error for values that aren't crazy tiny, since - // their relative error isn't representative. - if y.abs() > 0.0001 { - let x = i as f64 * lin_norm; - let e = (log_to_lin(x, offset, slope, params[0], params[1]) - y).abs() - / y.abs(); + for (x, y) in coords.iter().copied() { + // Only record error for the log part of the curve because we + // computed the linear segment's slope and offset analytically, + // and the relative error of points very near zero isn't reliable. + if y > transition.0 { + let e = (log_to_lin(x, offset, slope, log_offset, base) - y).abs() / y.abs(); max_err = max_err.max(e); avg_err += e; avg_samples += 1; @@ -47,11 +54,57 @@ pub fn find_parameters(lut: &[f32]) { } avg_err /= avg_samples as f64; - println!("Max Err: {:.4}%\nAvg Err: {:.4}%", max_err * 100.0, avg_err * 100.0); + println!( + "Max Err: {:.4}%\nAvg Err: {:.4}%", + max_err * 100.0, + avg_err * 100.0 + ); + + // dbg!(offset, log_offset, slope, base, end); println!( "{}{}", - crate::linear_log::generate_linear_to_log(offset, slope, params[0], params[1],), - crate::linear_log::generate_log_to_linear(offset, slope, params[0], params[1],), + crate::linear_log::generate_linear_to_log(offset, slope, log_offset, base), + crate::linear_log::generate_log_to_linear(offset, slope, log_offset, base), ); } + +/// This finds the minimum of functions with only one minimum (i.e. has +/// no local minimums other than the global one). It will not work +/// for functions that don't meet that criteria. +/// +/// It works by progressively narrowing the search range by: +/// 1. Splitting the range into four equal segments. +/// 2. Checking the slope of each segment. +/// 3. Narrowing the range to the two adjecent segments where +/// there is a switch from negative to positive slope. +fn optimize f64>(f: F, range: [f64; 2], iterations: usize) -> f64 { + let mut range = range; + + for _ in 0..iterations { + const SEG_POINTS: usize = 5; + let point = |xi| { + let n = xi as f64 / (SEG_POINTS - 1) as f64; + (range[0] * (1.0 - n)) + (range[1] * n) + }; + let mut last_xi = 0; + for i in 0..(SEG_POINTS - 1) { + last_xi = i; + let y1 = f(point(i)); + let y2 = f(point(i + 1)); + if (y2 - y1) >= 0.0 { + break; + } + } + let (r1, r2) = if last_xi == 0 { + (point(0), point(1)) + } else if last_xi == (SEG_POINTS - 1) { + (point(SEG_POINTS - 2), point(SEG_POINTS - 1)) + } else { + (point(last_xi - 1), point(last_xi + 1)) + }; + range = [r1, r2]; + } + + (range[0] + range[1]) * 0.5 +}