Merge pull request #16 from awxkee/dev

A bit changes

Author: awxkee

Cargo.toml

@@ -2,7 +2,7 @@ workspace = { members = ["app"] }
 
 [package]
 name = "gainforge"
-version = "0.3.7"
+version = "0.4.0"
 edition = "2021"
 description = "HDR tonemapping library"
 readme = "README.md"

app/src/main.rs

@@ -164,7 +164,7 @@ fn extract_images(file_path: &str) -> GainMapAssociationGroup {
 }
 
 fn main() {
-    let img = image::ImageReader::open("./assets/03.jpg")
+    let img = image::ImageReader::open("./assets/base.jpeg")
         .unwrap()
         .decode()
         .unwrap();
@@ -186,16 +186,16 @@ fn main() {
         //     saturation: Rgb::new(1.4, 1.4, 1.4),
         //     offset: Rgb::default(),
         // })),
-        ToneMappingMethod::Rec2408(GainHdrMetadata::new(2000f32, 203.)),
+        ToneMappingMethod::TunedReinhard(GainHdrMetadata::new(2000., 250.)),
         // ToneMappingMethod::ExtendedReinhard,
-        MappingColorSpace::Rgb(RgbToneMapperParameters {
-            gamut_clipping: GamutClipping::NoClip,
-            exposure: 1f32,
-        }),
-        // MappingColorSpace::YRgb(CommonToneMapperParameters {
-        //     exposure: 1f32,
+        // MappingColorSpace::Rgb(RgbToneMapperParameters {
         //     gamut_clipping: GamutClipping::NoClip,
+        //     exposure: 1f32,
         // }),
+        MappingColorSpace::Yrg(CommonToneMapperParameters {
+            exposure: 1f32,
+            gamut_clipping: GamutClipping::NoClip,
+        }),
         // MappingColorSpace::Jzazbz(JzazbzToneMapperParameters {
         //     content_brightness: 2000.,
         //     exposure: 1f32,
@@ -247,7 +247,7 @@ fn main() {
     // let compressed = dst.iter().map(|&x| (x >> 8) as u8).collect::<Vec<_>>();
 
     image::save_buffer(
-        "clamp_agx_t.jpg",
+        "tuned_reinhard.jpg",
         &dst,
         img.width(),
         img.height(),

src/err.rs

@@ -36,7 +36,7 @@ pub struct MismatchedSize {
     pub received: usize,
 }
 
-#[derive(Debug, Copy, Clone, Ord, PartialOrd, Eq, PartialEq)]
+#[derive(Debug, Clone, Ord, PartialOrd, Eq, PartialEq)]
 pub enum ForgeError {
     LaneSizeMismatch,
     LaneMultipleOfChannels,
@@ -52,6 +52,7 @@ pub enum ForgeError {
     InvalidGainMapConfiguration,
     ImageSizeMismatch,
     ZeroBaseSize,
+    UnsupportedToneMappingConfiguration(String),
     MinimumSliceSizeMismatch(MismatchedSize),
     MinimumStrideSizeMismatch(MismatchedSize),
     UnknownError,
@@ -90,6 +91,7 @@ impl Display for ForgeError {
                 "Minimum stride must have size at least {} but it is {}",
                 size.expected, size.received
             )),
+            ForgeError::UnsupportedToneMappingConfiguration(str) => f.write_str(str),
             ForgeError::UnknownError => f.write_str("Unknown error"),
         }
     }

src/gamma.rs

@@ -171,16 +171,27 @@ pub(crate) fn gamma2p8_to_linear(gamma: f32) -> f32 {
 
 #[inline(always)]
 /// Linear transfer function for PQ
-pub(crate) fn pq_to_linear(gamma: f32) -> f32 {
+pub(crate) fn pq_to_linear(gamma: f32, reference_display: f32) -> f32 {
     if gamma > 0.0 {
         let pow_gamma = f_powf(gamma, 1.0 / 78.84375);
         let num = (pow_gamma - 0.8359375).max(0.);
         let den = (18.8515625 - 18.6875 * pow_gamma).max(f32::MIN);
         let linear = f_powf(num / den, 1.0 / 0.1593017578125);
         // Scale so that SDR white is 1.0 (extended SDR).
-        const PQ_MAX_NITS: f32 = 10000.;
-        const SDR_WHITE_NITS: f32 = 203.;
-        linear * PQ_MAX_NITS / SDR_WHITE_NITS
+        linear * reference_display
+    } else {
+        0.0
+    }
+}
+
+#[inline(always)]
+/// Linear transfer function for PQ
+pub(crate) fn pq_to_linear_unscaled(gamma: f32) -> f32 {
+    if gamma > 0.0 {
+        let pow_gamma = f_powf(gamma, 1.0 / 78.84375);
+        let num = (pow_gamma - 0.8359375).max(0.);
+        let den = (18.8515625 - 18.6875 * pow_gamma).max(f32::MIN);
+        f_powf(num / den, 1.0 / 0.1593017578125)
     } else {
         0.0
     }
@@ -195,7 +206,22 @@ const HLG_WHITE_NITS: f32 = 1000.;
 pub(crate) fn pq_from_linear(linear: f32) -> f32 {
     if linear > 0.0 {
         // Scale from extended SDR range to [0.0, 1.0].
-        let linear = (linear * SDR_REFERENCE_DISPLAY / PQ_MAX_NITS).clamp(0., 1.);
+        let linear = (linear * (SDR_REFERENCE_DISPLAY / PQ_MAX_NITS)).clamp(0., 1.);
+        let pow_linear = f_powf(linear, 0.1593017578125);
+        let num = 0.1640625 * pow_linear - 0.1640625;
+        let den = 1.0 + 18.6875 * pow_linear;
+        f_powf(1.0 + num / den, 78.84375)
+    } else {
+        0.0
+    }
+}
+
+#[inline(always)]
+/// Gamma transfer function for PQ
+pub(crate) fn pq_from_linear_with_reference_display(linear: f32, reference_display: f32) -> f32 {
+    if linear > 0.0 {
+        // Scale from extended SDR range to [0.0, 1.0].
+        let linear = (linear * (reference_display * (1. / PQ_MAX_NITS))).clamp(0., 1.);
         let pow_linear = f_powf(linear, 0.1593017578125);
         let num = 0.1640625 * pow_linear - 0.1640625;
         let den = 1.0 + 18.6875 * pow_linear;
@@ -286,7 +312,7 @@ impl From<u8> for TransferFunction {
 
 impl TransferFunction {
     #[inline(always)]
-    pub fn linearize(&self, v: f32) -> f32 {
+    pub fn linearize(&self, v: f32, reference_display: f32) -> f32 {
         match self {
             TransferFunction::Srgb => srgb_to_linear(v),
             TransferFunction::Rec709 => rec709_to_linear(v),
@@ -296,7 +322,7 @@ impl TransferFunction {
             TransferFunction::Bt1361 => bt1361_to_linear(v),
             TransferFunction::Linear => trc_linear(v),
             TransferFunction::HybridLogGamma => hlg_to_linear(v),
-            TransferFunction::PerceptualQuantizer => pq_to_linear(v),
+            TransferFunction::PerceptualQuantizer => pq_to_linear(v, reference_display),
         }
     }
 
@@ -332,20 +358,24 @@ impl TransferFunction {
         table
     }
 
-    pub(crate) fn generate_linear_table_u16(&self, bit_depth: usize) -> Box<[f32; 65536]> {
+    pub(crate) fn generate_linear_table_u16(
+        &self,
+        bit_depth: usize,
+        reference_display: f32,
+    ) -> Box<[f32; 65536]> {
         let mut table = Box::new([0.; 65536]);
         let max_bp = (1 << bit_depth as u32) - 1;
         let max_scale = 1f32 / max_bp as f32;
         for (i, value) in table.iter_mut().take(max_bp).enumerate() {
-            *value = self.linearize(i as f32 * max_scale);
+            *value = self.linearize(i as f32 * max_scale, reference_display);
         }
         table
     }
 
-    pub(crate) fn generate_linear_table_u8(&self) -> Box<[f32; 256]> {
+    pub(crate) fn generate_linear_table_u8(&self, reference_display: f32) -> Box<[f32; 256]> {
         let mut table = Box::new([0.; 256]);
         for (i, value) in table.iter_mut().enumerate() {
-            *value = self.linearize(i as f32 / 255.);
+            *value = self.linearize(i as f32 * (1. / 255.), reference_display);
         }
         table
     }

src/lib.rs

@@ -257,7 +257,7 @@ pub use iso_gain_map::{
     UhdrDirectorySeq, UhdrItem, UhdrItemContainerLi, UhdrItemResource,
 };
 pub use mappers::{AgxCustomLook, AgxLook, ToneMappingMethod};
-use num_traits::Num;
+use num_traits::{Float, Num};
 pub use spline::FilmicSplineParameters;
 pub use tonemapper::{
     create_tone_mapper_rgb, create_tone_mapper_rgb10, create_tone_mapper_rgb12,
@@ -269,13 +269,6 @@ pub use tonemapper::{
 };
 
 #[inline]
-pub(crate) fn m_clamp<T: Num + PartialOrd>(a: T, min: T, max: T) -> T {
-    if a > max {
-        max
-    } else if a >= min {
-        a
-    } else {
-        // a < min or a is NaN
-        min
-    }
+pub(crate) fn m_clamp<T: Num + Float>(a: T, min: T, max: T) -> T {
+    a.min(max).max(min)
 }

src/mappers.rs

@@ -26,6 +26,7 @@
  * // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
  * // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  */
+use crate::gamma::{pq_from_linear_with_reference_display, pq_to_linear_unscaled};
 use crate::mlaf::{fmla, mlaf};
 use crate::spline::FilmicSplineParameters;
 use crate::GainHdrMetadata;
@@ -49,8 +50,10 @@ pub enum AgxLook {
 /// many of the supported tone mapping methods.
 #[derive(Debug, Copy, Clone, PartialOrd, PartialEq)]
 pub enum ToneMappingMethod {
+    /// Fast and accurate tuned Reinhard, preferred
+    TunedReinhard(GainHdrMetadata),
     /// ITU-R broadcasting TV [recommendation 2408](https://www.itu.int/dms_pub/itu-r/opb/rep/R-REP-BT.2408-4-2021-PDF-E.pdf)
-    Rec2408(GainHdrMetadata),
+    Itu2408(GainHdrMetadata),
     /// The ['Uncharted 2' filmic](https://www.gdcvault.com/play/1012351/Uncharted-2-HDR)
     /// tone mapping method.
     Filmic,
@@ -84,12 +87,146 @@ pub(crate) trait ToneMap {
 
 #[derive(Debug, Clone, Copy)]
 pub(crate) struct Rec2408ToneMapper<const CN: usize> {
+    content_max_brightness: f32,
+    display_max_brightness: f32,
+    primaries: [f32; 3],
+    content_min_luminance: f32,
+    content_luminance_range: f32,
+    inv_pq_mastering_range: f32,
+    min_lum: f32,
+    max_lum: f32,
+    ks: f32,
+    inv_one_minus_ks: f32,
+    one_minus_ks: f32,
+    normalizer: f32,
+    inv_target_peak: f32,
+}
+
+impl<const CN: usize> Rec2408ToneMapper<CN> {
+    pub(crate) fn new(
+        content_max_brightness: f32,
+        display_max_brightness: f32,
+        primaries: [f32; 3],
+    ) -> Self {
+        let content_luminance_black = pq_from_linear_with_reference_display(0., 10000.0);
+        let content_luminance_white =
+            pq_from_linear_with_reference_display(content_max_brightness / 10000.0, 10000.0);
+        let content_luminance_range = content_luminance_white - content_luminance_black;
+        let inv_pq_mastering_range = 1.0 / content_luminance_range;
+
+        let min_lum = (pq_from_linear_with_reference_display(0., 10000.0)
+            - content_luminance_black)
+            * inv_pq_mastering_range;
+        let max_lum =
+            (pq_from_linear_with_reference_display(display_max_brightness / 10000.0, 10000.0)
+                - content_luminance_black)
+                * inv_pq_mastering_range;
+        let ks = 1.5 * max_lum - 0.5;
+
+        Self {
+            content_max_brightness,
+            display_max_brightness,
+            primaries,
+            content_min_luminance: content_luminance_black,
+            content_luminance_range,
+            inv_pq_mastering_range,
+            min_lum,
+            max_lum,
+            ks,
+            inv_one_minus_ks: 1.0 / (1.0 - ks).max(1e-6),
+            one_minus_ks: 1.0 - ks,
+            normalizer: content_max_brightness / display_max_brightness,
+            inv_target_peak: 1.0 / display_max_brightness,
+        }
+    }
+}
+
+impl<const CN: usize> Rec2408ToneMapper<CN> {
+    #[inline(always)]
+    fn t(&self, a: f32) -> f32 {
+        (a - self.ks) * self.inv_one_minus_ks
+    }
+
+    #[inline]
+    fn hermite_spline(&self, b: f32) -> f32 {
+        let t_b = self.t(b);
+        let t_b_2 = t_b * t_b;
+        let t_b_3 = t_b_2 * t_b;
+        fmla(2.0, t_b_3, fmla(-3.0, t_b_2, 1.0)) * self.ks
+            + fmla(-2.0, t_b_2, t_b_3 + t_b) * self.one_minus_ks
+            + fmla(-2.0, t_b_3, 3.0 * t_b_2) * self.max_lum
+    }
+
+    #[inline(always)]
+    fn make_luma_scale(&self, luma: f32) -> f32 {
+        let s = pq_from_linear_with_reference_display(luma / 10000.0, 10000.0);
+        let normalized_pq =
+            ((s - self.content_min_luminance) * self.inv_pq_mastering_range).min(1.0);
+
+        let e2 = if normalized_pq < self.ks {
+            normalized_pq
+        } else {
+            self.hermite_spline(normalized_pq)
+        };
+
+        let one_minus_e2 = 1.0 - e2;
+        let one_minus_e2_2 = one_minus_e2 * one_minus_e2;
+        let one_minus_e2_4 = one_minus_e2_2 * one_minus_e2_2;
+        let e3 = fmla(self.min_lum, one_minus_e2_4, e2);
+        let e4 = e3 * self.content_luminance_range + self.content_min_luminance;
+        let d4 = pq_to_linear_unscaled(e4) * 10000.0;
+        let new_luminance = d4.min(self.display_max_brightness).max(0.);
+
+        let min_luminance = 1e-6;
+        let use_limit = luma <= min_luminance;
+        let ratio = new_luminance / luma.max(min_luminance);
+        let limit = new_luminance * self.inv_target_peak;
+        let scale = ratio * self.normalizer;
+        if use_limit {
+            limit
+        } else {
+            scale
+        }
+    }
+}
+
+impl<const CN: usize> ToneMap for Rec2408ToneMapper<CN> {
+    fn process_lane(&self, in_place: &mut [f32]) {
+        for chunk in in_place.chunks_exact_mut(CN) {
+            let luma = fmla(
+                chunk[0],
+                self.primaries[0],
+                fmla(chunk[1], self.primaries[1], chunk[2] * self.primaries[2]),
+            ) * self.content_max_brightness;
+            if luma == 0. {
+                chunk[0] = 0.;
+                chunk[1] = 0.;
+                chunk[2] = 0.;
+                continue;
+            }
+            let scale = self.make_luma_scale(luma);
+            chunk[0] *= scale;
+            chunk[1] *= scale;
+            chunk[2] *= scale;
+        }
+    }
+
+    fn process_luma_lane(&self, chunk: &mut [f32]) {
+        for chunk in chunk.chunks_exact_mut(CN) {
+            let scale = self.make_luma_scale(chunk[0] * self.content_max_brightness);
+            chunk[0] *= scale;
+        }
+    }
+}
+
+#[derive(Debug, Clone, Copy)]
+pub(crate) struct TunedReinhardToneMapper<const CN: usize> {
     w_a: f32,
     w_b: f32,
     primaries: [f32; 3],
 }
 
-impl<const CN: usize> Rec2408ToneMapper<CN> {
+impl<const CN: usize> TunedReinhardToneMapper<CN> {
     pub(crate) fn new(
         content_max_brightness: f32,
         display_max_brightness: f32,
@@ -107,14 +244,14 @@ impl<const CN: usize> Rec2408ToneMapper<CN> {
     }
 }
 
-impl<const CN: usize> Rec2408ToneMapper<CN> {
+impl<const CN: usize> TunedReinhardToneMapper<CN> {
     #[inline(always)]
     fn tonemap(&self, luma: f32) -> f32 {
         mlaf(1f32, self.w_a, luma) / mlaf(1f32, self.w_b, luma)
     }
 }
 
-impl<const CN: usize> ToneMap for Rec2408ToneMapper<CN> {
+impl<const CN: usize> ToneMap for TunedReinhardToneMapper<CN> {
     fn process_lane(&self, in_place: &mut [f32]) {
         for chunk in in_place.chunks_exact_mut(CN) {
             let luma = fmla(