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

@@ -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"

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"

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

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<f32> = 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()),

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.
+    // Collect LUT points as (x, y) coordinates.
-    let idxs: Vec<_> = (0..lut.len()).step_by(lut.len() / 256).collect();
+    let coords: Vec<(f64, f64)> = lut
-    let coords: Vec<(f64, f64)> = idxs
         .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 base = optimize(
-        let mut avg_sqr_err = 0.0f64;
+        |v: f64| {
-        for (x, y) in coords.iter().copied() {
+            let log_offset = crate::linear_log::find_log_offset_for_end(end, offset, slope, v);
-            let e = (log_to_lin(x, offset, slope, v[0], v[1]) - y).abs() / y.abs();
+            let mut sqr_err = 0.0f64;
-            avg_sqr_err += e * e;
+            for (x, y) in coords.iter().copied() {
-        }
+                let e = (log_to_lin(x, offset, slope, log_offset, v) - y).abs() / y.abs();
-        let last_y = lut[lut.len() - 1] as f64;
+                sqr_err += e * e;
-        let e = (log_to_lin(1.0, offset, slope, v[0], v[1]) - last_y).abs() / last_y.abs();
+            }
-        avg_sqr_err += e * e;
+            sqr_err
-        avg_sqr_err
+        },
-    };
+        [1.5, 10000000.0],
-    let input_interval = vec![(-0.2, 0.2), (1.1, 1000.0)];
+        100,
-    let (_, params) = simplers_optimization::Optimizer::minimize(&f, &input_interval, 1000000);
+    );
+    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() {
+    for (x, y) in coords.iter().copied() {
-        // We only record error for values that aren't crazy tiny, since
+        // Only record error for the log part of the curve because we
-        // their relative error isn't representative.
+        // computed the linear segment's slope and offset analytically,
-        if y.abs() > 0.0001 {
+        // and the relative error of points very near zero isn't reliable.
-            let x = i as f64 * lin_norm;
+        if y > transition.0 {
-            let e = (log_to_lin(x, offset, slope, params[0], params[1]) - y).abs()
+            let e = (log_to_lin(x, offset, slope, log_offset, base) - y).abs() / 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_linear_to_log(offset, slope, log_offset, base),
-        crate::linear_log::generate_log_to_linear(offset, slope, params[0], params[1],),
+        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<F: Fn(f64) -> 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
+}