gainforge 0.5.0

/*
 * // Copyright (c) Radzivon Bartoshyk 2/2025. All rights reserved.
 * //
 * // Redistribution and use in source and binary forms, with or without modification,
 * // are permitted provided that the following conditions are met:
 * //
 * // 1.  Redistributions of source code must retain the above copyright notice, this
 * // list of conditions and the following disclaimer.
 * //
 * // 2.  Redistributions in binary form must reproduce the above copyright notice,
 * // this list of conditions and the following disclaimer in the documentation
 * // and/or other materials provided with the distribution.
 * //
 * // 3.  Neither the name of the copyright holder nor the names of its
 * // contributors may be used to endorse or promote products derived from
 * // this software without specific prior written permission.
 * //
 * // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
 * // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
 * // 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::ToneMappingMethod;
use crate::err::ForgeError;
use crate::gamma::trc_from_cicp;
use crate::mappers::{
    AcesToneMapper, AgxDefault, AgxGolden, AgxLook, AgxPunchy, AgxToneMapper, ClampToneMapper,
    FilmicToneMapper, Rec2408ToneMapper, ReinhardJodieToneMapper, ReinhardToneMapper, ToneMap,
};
use crate::mlaf::fmla;
use crate::rgb_tone_mapper::{MatrixGamutClipping, MatrixStage, ToneMapperImpl};
use crate::spline::{SplineToneMapper, create_spline};
#[allow(unused_imports)]
use crate::util::MangledCoercion;
use moxcms::{ColorProfile, InPlaceStage, Jzazbz, Matrix3f, Rgb, Yrg, filmlike_clip};
use num_traits::AsPrimitive;
use std::fmt::Debug;
use std::sync::Arc;

/// Defines gamut clipping mode
#[derive(Debug, Default, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, Hash)]
pub enum GamutClipping {
    #[default]
    NoClip,
    Clip,
}

#[derive(Debug, Copy, Clone, PartialOrd, PartialEq)]
pub enum MappingColorSpace {
    /// Scene linear light RGB colorspace.
    ///
    /// Some description what to expect:
    /// - Does not produce colours outside the destination gamut.
    /// - Generally tends to produce “natural” looking colours, although saturation of extreme colours
    ///   is reduced substantially and some colours may be changed in hue.
    /// - Problematic when it is desired to retain bright saturated colours, such as coloured lights at
    ///   night
    /// - Depending on the amount of compression, the saturation decrease may be excessive, and
    ///   occasionally hue changes can be objectionable.
    Rgb(RgbToneMapperParameters),
    /// Yrg filmic colorspace.
    ///
    /// Slightly better results as in linear RGB, for a little computational cost.
    ///
    /// Some description what to expect:
    /// - Has the potential to produce colours outside the destination gamut, which then require gamut
    ///   mapping
    /// - Preserves chromaticity except for where gamut mapping is applied. Does not produce a
    ///   “natural” looking desaturation of tonally compressed colours.
    /// - Problems can be avoided by using in combination with a variable desaturation and gamut
    ///   mapping algorithm, although such algorithms generally perform best in hue, saturation and
    ///   lightness colour spaces (requiring a colour space change).
    ///
    /// Algorithm here do not perform gamut mapping.
    Yrg(CommonToneMapperParameters),
    /// JzAzBz perceptual HDR colorspace.
    ///
    /// Some description of what to expect:
    /// - Designed for HDR workflows, explicitly modeling display-referred luminance.
    /// - Thus, as at the first point often produces the best results.
    /// - Provides better perceptual uniformity than Oklab, particularly in high dynamic range content.
    /// - Preserves hue well, but chroma adjustments may be necessary when mapping between displays with different peak brightness.
    /// - Can help avoid color distortions and oversaturation when adapting HDR content to lower-luminance displays.
    ///
    /// This is very slow computational method.
    /// Content brightness should be specified.
    Jzazbz(JzazbzToneMapperParameters),
}

#[derive(Debug, Copy, Clone, PartialOrd, PartialEq)]
pub struct JzazbzToneMapperParameters {
    pub content_brightness: f32,
    pub exposure: f32,
    pub gamut_clipping: GamutClipping,
}

#[derive(Debug, Copy, Clone, PartialOrd, PartialEq)]
pub struct CommonToneMapperParameters {
    pub exposure: f32,
    pub gamut_clipping: GamutClipping,
}

#[derive(Debug, Copy, Clone, PartialOrd, PartialEq)]
pub struct RgbToneMapperParameters {
    pub exposure: f32,
    pub gamut_clipping: GamutClipping,
}

impl Default for JzazbzToneMapperParameters {
    fn default() -> Self {
        Self {
            content_brightness: 1000f32,
            exposure: 1.0f32,
            gamut_clipping: GamutClipping::Clip,
        }
    }
}

impl Default for CommonToneMapperParameters {
    fn default() -> Self {
        Self {
            exposure: 1.0f32,
            gamut_clipping: GamutClipping::Clip,
        }
    }
}

impl Default for RgbToneMapperParameters {
    fn default() -> Self {
        Self {
            exposure: 1.0f32,
            gamut_clipping: GamutClipping::default(),
        }
    }
}

pub type ToneMapping<T> = Arc<dyn ToneMapper<T> + Send + Sync>;

pub(crate) type SyncToneMap = dyn ToneMap + Send + Sync;

struct ToneMapperImplYrg<T: Copy, const N: usize, const CN: usize, const GAMMA_SIZE: usize> {
    linear_map_r: Box<[f32; N]>,
    linear_map_g: Box<[f32; N]>,
    linear_map_b: Box<[f32; N]>,
    gamma_map_r: Box<[T; 65536]>,
    gamma_map_g: Box<[T; 65536]>,
    gamma_map_b: Box<[T; 65536]>,
    to_xyz: Matrix3f,
    to_rgb: Matrix3f,
    tone_map: Arc<SyncToneMap>,
    parameters: CommonToneMapperParameters,
}

struct ToneMapperImplJzazbz<T: Copy, const N: usize, const CN: usize, const GAMMA_SIZE: usize> {
    linear_map_r: Box<[f32; N]>,
    linear_map_g: Box<[f32; N]>,
    linear_map_b: Box<[f32; N]>,
    gamma_map_r: Box<[T; 65536]>,
    gamma_map_g: Box<[T; 65536]>,
    gamma_map_b: Box<[T; 65536]>,
    to_xyz: Matrix3f,
    to_rgb: Matrix3f,
    tone_map: Arc<SyncToneMap>,
    parameters: JzazbzToneMapperParameters,
}

pub trait ToneMapper<T: Copy + Default + Debug> {
    /// Tone map image lane.
    ///
    /// Lane length must be multiple of channels.
    /// Lane length must match.
    fn tonemap_lane(&self, src: &[T], dst: &mut [T]) -> Result<(), ForgeError>;

    /// Tone map lane whereas content been linearized.
    ///
    /// Lane length must be multiple of channels.
    fn tonemap_linearized_lane(&self, in_place: &mut [f32]) -> Result<(), ForgeError>;
}

impl<
    T: Copy + AsPrimitive<usize> + Clone + Default + Debug,
    const N: usize,
    const CN: usize,
    const GAMMA_SIZE: usize,
> ToneMapper<T> for ToneMapperImplYrg<T, N, CN, GAMMA_SIZE>
where
    u32: AsPrimitive<T>,
{
    fn tonemap_lane(&self, src: &[T], dst: &mut [T]) -> Result<(), ForgeError> {
        assert!(CN == 3 || CN == 4);
        if src.len() != dst.len() {
            return Err(ForgeError::LaneSizeMismatch);
        }
        if !src.len().is_multiple_of(CN) {
            return Err(ForgeError::LaneMultipleOfChannels);
        }
        assert_eq!(src.len(), dst.len());
        let mut linearized_content = vec![0f32; src.len()];
        for (src, dst) in src
            .as_chunks::<CN>()
            .0
            .iter()
            .zip(linearized_content.as_chunks_mut::<CN>().0.iter_mut())
        {
            let xyz = (Rgb::new(
                self.linear_map_r[src[0].as_()],
                self.linear_map_g[src[1].as_()],
                self.linear_map_b[src[2].as_()],
            ) * self.parameters.exposure)
                .to_xyz(self.to_xyz);
            let yrg = Yrg::from_xyz(xyz);
            dst[0] = yrg.y;
            dst[1] = yrg.r;
            dst[2] = yrg.g;
            if CN == 4 {
                dst[3] = f32::from_bits(src[3].as_() as u32);
            }
        }

        self.tonemap_linearized_lane(&mut linearized_content)?;

        match self.parameters.gamut_clipping {
            GamutClipping::NoClip => {
                for dst in linearized_content.as_chunks_mut::<CN>().0.iter_mut() {
                    let yrg = Yrg::new(dst[0], dst[1], dst[2]);
                    let xyz = yrg.to_xyz();
                    let rgb = xyz.to_linear_rgb(self.to_rgb);
                    dst[0] = rgb.r.min(1.).max(0.);
                    dst[1] = rgb.g.min(1.).max(0.);
                    dst[2] = rgb.b.min(1.).max(0.);
                }
            }
            GamutClipping::Clip => {
                for dst in linearized_content.as_chunks_mut::<CN>().0.iter_mut() {
                    let yrg = Yrg::new(dst[0], dst[1], dst[2]);
                    let xyz = yrg.to_xyz();
                    let mut rgb = xyz.to_linear_rgb(self.to_rgb);
                    if rgb.is_out_of_gamut() {
                        rgb = filmlike_clip(rgb);
                    }
                    dst[0] = rgb.r.min(1.).max(0.);
                    dst[1] = rgb.g.min(1.).max(0.);
                    dst[2] = rgb.b.min(1.).max(0.);
                }
            }
        }

        let scale_value = (GAMMA_SIZE - 1) as f32;

        for (dst, src) in dst
            .as_chunks_mut::<CN>()
            .0
            .iter_mut()
            .zip(linearized_content.as_chunks::<CN>().0.iter())
        {
            let r = fmla(src[0], scale_value, 0.5) as u16;
            let g = fmla(src[1], scale_value, 0.5) as u16;
            let b = fmla(src[2], scale_value, 0.5) as u16;
            dst[0] = self.gamma_map_r[r as usize];
            dst[1] = self.gamma_map_g[g as usize];
            dst[2] = self.gamma_map_b[b as usize];
            if CN == 4 {
                dst[3] = src[3].to_bits().as_();
            }
        }

        Ok(())
    }

    fn tonemap_linearized_lane(&self, in_place: &mut [f32]) -> Result<(), ForgeError> {
        assert!(CN == 3 || CN == 4);
        if !in_place.len().is_multiple_of(CN) {
            return Err(ForgeError::LaneMultipleOfChannels);
        }
        self.tone_map.process_luma_lane(in_place);
        Ok(())
    }
}

impl<
    T: Copy + AsPrimitive<usize> + Clone + Default + Debug,
    const N: usize,
    const CN: usize,
    const GAMMA_SIZE: usize,
> ToneMapper<T> for ToneMapperImplJzazbz<T, N, CN, GAMMA_SIZE>
where
    u32: AsPrimitive<T>,
{
    fn tonemap_lane(&self, src: &[T], dst: &mut [T]) -> Result<(), ForgeError> {
        assert!(CN == 3 || CN == 4);
        if src.len() != dst.len() {
            return Err(ForgeError::LaneSizeMismatch);
        }
        if !src.len().is_multiple_of(CN) {
            return Err(ForgeError::LaneMultipleOfChannels);
        }
        assert_eq!(src.len(), dst.len());
        let mut linearized_content = vec![0f32; src.len()];

        for (src, dst) in src
            .as_chunks::<CN>()
            .0
            .iter()
            .zip(linearized_content.as_chunks_mut::<CN>().0.iter_mut())
        {
            let xyz = (Rgb::new(
                self.linear_map_r[src[0].as_()],
                self.linear_map_g[src[1].as_()],
                self.linear_map_b[src[2].as_()],
            ) * self.parameters.exposure)
                .to_xyz(self.to_xyz);
            let jab =
                Jzazbz::from_xyz_with_display_luminance(xyz, self.parameters.content_brightness);
            dst[0] = jab.jz;
            dst[1] = jab.az;
            dst[2] = jab.bz;
            if CN == 4 {
                dst[3] = f32::from_bits(src[3].as_() as u32);
            }
        }

        self.tonemap_linearized_lane(&mut linearized_content)?;

        match self.parameters.gamut_clipping {
            GamutClipping::NoClip => {
                for dst in linearized_content.as_chunks_mut::<CN>().0.iter_mut() {
                    let jab = Jzazbz::new(dst[0], dst[1], dst[2]);
                    let xyz = jab.to_xyz(self.parameters.content_brightness);
                    let rgb = xyz.to_linear_rgb(self.to_rgb);
                    dst[0] = rgb.r.min(1.).max(0.);
                    dst[1] = rgb.g.min(1.).max(0.);
                    dst[2] = rgb.b.min(1.).max(0.);
                }
            }
            GamutClipping::Clip => {
                for dst in linearized_content.as_chunks_mut::<CN>().0.iter_mut() {
                    let jab = Jzazbz::new(dst[0], dst[1], dst[2]);
                    let xyz = jab.to_xyz(self.parameters.content_brightness);
                    let mut rgb = xyz.to_linear_rgb(self.to_rgb);
                    if rgb.is_out_of_gamut() {
                        rgb = filmlike_clip(rgb);
                    }
                    dst[0] = rgb.r.min(1.).max(0.);
                    dst[1] = rgb.g.min(1.).max(0.);
                    dst[2] = rgb.b.min(1.).max(0.);
                }
            }
        }

        let scale_value = (GAMMA_SIZE - 1) as f32;

        for (dst, src) in dst
            .as_chunks_mut::<CN>()
            .0
            .iter_mut()
            .zip(linearized_content.as_chunks::<CN>().0.iter())
        {
            let r = fmla(src[0], scale_value, 0.5) as u16;
            let g = fmla(src[1], scale_value, 0.5) as u16;
            let b = fmla(src[2], scale_value, 0.5) as u16;
            dst[0] = self.gamma_map_r[r as usize];
            dst[1] = self.gamma_map_g[g as usize];
            dst[2] = self.gamma_map_b[b as usize];
            if CN == 4 {
                dst[3] = src[3].to_bits().as_();
            }
        }

        Ok(())
    }

    fn tonemap_linearized_lane(&self, in_place: &mut [f32]) -> Result<(), ForgeError> {
        assert!(CN == 3 || CN == 4);
        if !in_place.len().is_multiple_of(CN) {
            return Err(ForgeError::LaneMultipleOfChannels);
        }
        self.tone_map.process_luma_lane(in_place);
        Ok(())
    }
}

#[derive(Debug, Clone, Copy, PartialOrd, PartialEq)]
pub struct GainHdrMetadata {
    pub content_max_brightness: f32,
    pub display_max_brightness: f32,
}

impl Default for GainHdrMetadata {
    fn default() -> Self {
        Self {
            content_max_brightness: 1000f32,
            display_max_brightness: 250f32,
        }
    }
}

impl GainHdrMetadata {
    pub fn new(content_max_brightness: f32, display_max_brightness: f32) -> Self {
        Self {
            content_max_brightness,
            display_max_brightness,
        }
    }
}

fn make_icc_transform(
    input_color_space: &ColorProfile,
    output_color_space: &ColorProfile,
) -> Matrix3f {
    input_color_space
        .transform_matrix(output_color_space)
        .to_f32()
}

fn matrix_stage<const CN: usize>(
    gamut_clipping: GamutClipping,
    conversion: Matrix3f,
) -> Box<dyn InPlaceStage + Send + Sync> {
    if gamut_clipping == GamutClipping::Clip {
        Box::new(MatrixGamutClipping::<CN> {
            gamut_color_conversion: conversion,
        })
    } else {
        #[cfg(target_arch = "x86_64")]
        {
            if std::arch::is_x86_feature_detected!("avx2")
                && std::arch::is_x86_feature_detected!("fma")
            {
                use crate::rgb_tone_mapper::FmaMatrixStage;
                return Box::new(FmaMatrixStage::<CN> {
                    gamut_color_conversion: conversion,
                });
            }
        }
        Box::new(MatrixStage::<CN> {
            gamut_color_conversion: conversion,
        })
    }
}

fn create_tone_mapper_u8<const CN: usize>(
    input_color_space: &ColorProfile,
    output_color_space: &ColorProfile,
    method: ToneMappingMethod,
    working_color_space: MappingColorSpace,
) -> Result<ToneMapping<u8>, ForgeError> {
    let (linear_table_r, linear_table_g, linear_table_b);
    if let Some(trc) = input_color_space
        .cicp
        .and_then(|x| trc_from_cicp(x.transfer_characteristics))
    {
        const PQ_MAX_NITS: f32 = 10000.;
        const SDR_WHITE_NITS: f32 = 203.;
        let reference_display = PQ_MAX_NITS / SDR_WHITE_NITS;
        linear_table_r =
            trc.generate_linear_table_u8(if matches!(method, ToneMappingMethod::Itu2408(_)) {
                1.
            } else {
                reference_display
            });
        linear_table_g = linear_table_r.clone();
        linear_table_b = linear_table_g.clone();
    } else {
        linear_table_r = input_color_space
            .build_r_linearize_table::<u8, 256, 8>(true)
            .map_err(|_| ForgeError::InvalidTrcCurve)?;
        linear_table_g = input_color_space
            .build_g_linearize_table::<u8, 256, 8>(true)
            .map_err(|_| ForgeError::InvalidTrcCurve)?;
        linear_table_b = input_color_space
            .build_b_linearize_table::<u8, 256, 8>(true)
            .map_err(|_| ForgeError::InvalidTrcCurve)?;
    }
    let (gamma_table_r, gamma_table_g, gamma_table_b);
    if let Some(trc) = output_color_space
        .cicp
        .and_then(|x| trc_from_cicp(x.transfer_characteristics))
    {
        gamma_table_r = trc.generate_gamma_table_u8();
        gamma_table_g = gamma_table_r.clone();
        gamma_table_b = gamma_table_g.clone();
    } else {
        gamma_table_r = output_color_space
            .build_gamma_table::<u8, 65536, 4096, 8>(&output_color_space.red_trc, true)
            .unwrap();
        gamma_table_g = output_color_space
            .build_gamma_table::<u8, 65536, 4096, 8>(&output_color_space.green_trc, true)
            .unwrap();
        gamma_table_b = output_color_space
            .build_gamma_table::<u8, 65536, 4096, 8>(&output_color_space.blue_trc, true)
            .unwrap();
    }
    let conversion = make_icc_transform(input_color_space, output_color_space);

    let tone_map = make_mapper::<CN>(input_color_space, method);

    let _all_luts_the_same = linear_table_r == linear_table_g
        && linear_table_b == linear_table_g
        && gamma_table_r == gamma_table_g
        && gamma_table_g == gamma_table_b;

    match working_color_space {
        MappingColorSpace::Rgb(params) => {
            #[cfg(all(target_arch = "x86_64", feature = "avx"))]
            {
                if _all_luts_the_same
                    && let Some(value) = gen_avx_transform::<u8, 256, CN>(
                        input_color_space,
                        &method,
                        &working_color_space,
                        &linear_table_r,
                        &gamma_table_r,
                        &tone_map,
                        conversion,
                        8,
                        4096,
                    )
                {
                    return value;
                }
            }
            #[cfg(all(target_arch = "aarch64", feature = "neon"))]
            {
                if _all_luts_the_same
                    && let Some(value) = get_neon_transform::<u8, 256, CN>(
                        input_color_space,
                        &method,
                        &working_color_space,
                        &linear_table_r,
                        &gamma_table_r,
                        &tone_map,
                        conversion,
                        8,
                        4096,
                    )
                {
                    return value;
                }
            }
            let im_stage: Box<dyn InPlaceStage + Send + Sync> =
                matrix_stage::<CN>(params.gamut_clipping, conversion);

            Ok(Arc::new(ToneMapperImpl::<u8, 256, CN, 4096> {
                linear_map_r: linear_table_r,
                linear_map_g: linear_table_g,
                linear_map_b: linear_table_b,
                gamma_map_r: gamma_table_r,
                gamma_map_b: gamma_table_g,
                gamma_map_g: gamma_table_b,
                im_stage: Some(im_stage),
                tone_map,
                params,
            }))
        }
        MappingColorSpace::Yrg(params) => {
            // We need to adapt PCS D50 to CIE XYZ 2006 with D65 white point first.
            let mut to_xyz = input_color_space.rgb_to_xyz_matrix().to_f32();
            to_xyz = to_xyz.mul_row::<1>(1.05785528f32);
            let output_d50 = output_color_space.rgb_to_xyz_matrix().to_f32();
            let mut to_rgb = output_d50.inverse();
            to_rgb = to_rgb.mul_row::<1>(1. / 1.05785528f32);

            Ok(Arc::new(ToneMapperImplYrg::<u8, 256, CN, 4096> {
                linear_map_r: linear_table_r,
                linear_map_g: linear_table_g,
                linear_map_b: linear_table_b,
                gamma_map_r: gamma_table_r,
                gamma_map_b: gamma_table_g,
                gamma_map_g: gamma_table_b,
                to_xyz,
                to_rgb,
                tone_map,
                parameters: params,
            }))
        }
        MappingColorSpace::Jzazbz(brightness) => {
            // We need to adapt PCS D50 to XYZ with D65 white point first.
            let to_xyz = input_color_space.rgb_to_xyz_matrix().to_f32();
            let output_d65 = output_color_space.rgb_to_xyz_matrix().to_f32();
            let to_rgb = output_d65.inverse();

            Ok(Arc::new(ToneMapperImplJzazbz::<u8, 256, CN, 4096> {
                linear_map_r: linear_table_r,
                linear_map_g: linear_table_g,
                linear_map_b: linear_table_b,
                gamma_map_r: gamma_table_r,
                gamma_map_b: gamma_table_g,
                gamma_map_g: gamma_table_b,
                to_xyz,
                to_rgb,
                tone_map,
                parameters: brightness,
            }))
        }
    }
}

fn create_tone_mapper_u16<const CN: usize, const BIT_DEPTH: usize>(
    input_color_space: &ColorProfile,
    output_color_space: &ColorProfile,
    method: ToneMappingMethod,
    working_color_space: MappingColorSpace,
) -> Result<ToneMapping<u16>, ForgeError> {
    assert!((8..=16).contains(&BIT_DEPTH));
    let (linear_table_r, linear_table_g, linear_table_b);
    if let Some(trc) = input_color_space
        .cicp
        .and_then(|x| trc_from_cicp(x.transfer_characteristics))
    {
        const PQ_MAX_NITS: f32 = 10000.;
        const SDR_WHITE_NITS: f32 = 203.;
        let reference_display = PQ_MAX_NITS / SDR_WHITE_NITS;
        linear_table_r = trc.generate_linear_table_u16(
            BIT_DEPTH,
            if matches!(method, ToneMappingMethod::Itu2408(_)) {
                1.
            } else {
                reference_display
            },
        );
        linear_table_g = linear_table_r.clone();
        linear_table_b = linear_table_g.clone();
    } else {
        linear_table_r = input_color_space
            .build_r_linearize_table::<u16, 65536, BIT_DEPTH>(true)
            .map_err(|_| ForgeError::InvalidTrcCurve)?;
        linear_table_g = input_color_space
            .build_g_linearize_table::<u16, 65536, BIT_DEPTH>(true)
            .map_err(|_| ForgeError::InvalidTrcCurve)?;
        linear_table_b = input_color_space
            .build_b_linearize_table::<u16, 65536, BIT_DEPTH>(true)
            .map_err(|_| ForgeError::InvalidTrcCurve)?;
    }
    let (gamma_table_r, gamma_table_g, gamma_table_b);
    if let Some(trc) = output_color_space
        .cicp
        .and_then(|x| trc_from_cicp(x.transfer_characteristics))
    {
        gamma_table_r = trc.generate_gamma_table_u16(BIT_DEPTH);
        gamma_table_g = gamma_table_r.clone();
        gamma_table_b = gamma_table_g.clone();
    } else {
        gamma_table_r = output_color_space
            .build_gamma_table::<u16, 65536, 65536, BIT_DEPTH>(&output_color_space.red_trc, true)
            .unwrap();
        gamma_table_g = output_color_space
            .build_gamma_table::<u16, 65536, 65536, BIT_DEPTH>(&output_color_space.green_trc, true)
            .unwrap();
        gamma_table_b = output_color_space
            .build_gamma_table::<u16, 65536, 65536, BIT_DEPTH>(&output_color_space.blue_trc, true)
            .unwrap();
    }
    let tone_map = make_mapper::<CN>(input_color_space, method);

    let _all_luts_the_same = linear_table_r == linear_table_g
        && linear_table_b == linear_table_g
        && gamma_table_r == gamma_table_g
        && gamma_table_g == gamma_table_b;

    match working_color_space {
        MappingColorSpace::Rgb(params) => {
            let conversion = make_icc_transform(input_color_space, output_color_space);
            #[cfg(all(target_arch = "x86_64", feature = "avx"))]
            {
                if _all_luts_the_same
                    && let Some(value) = gen_avx_transform::<u16, 65536, CN>(
                        input_color_space,
                        &method,
                        &working_color_space,
                        &linear_table_r,
                        &gamma_table_r,
                        &tone_map,
                        conversion,
                        BIT_DEPTH,
                        65536,
                    )
                {
                    return value;
                }
            }
            #[cfg(all(target_arch = "aarch64", feature = "neon"))]
            {
                if _all_luts_the_same
                    && let Some(value) = get_neon_transform::<u16, 65536, CN>(
                        input_color_space,
                        &method,
                        &working_color_space,
                        &linear_table_r,
                        &gamma_table_r,
                        &tone_map,
                        conversion,
                        BIT_DEPTH,
                        65536,
                    )
                {
                    return value;
                }
            }

            let im_stage: Box<dyn InPlaceStage + Send + Sync> =
                matrix_stage::<CN>(params.gamut_clipping, conversion);

            Ok(Arc::new(ToneMapperImpl::<u16, 65536, CN, 65536> {
                linear_map_r: linear_table_r,
                linear_map_g: linear_table_g,
                linear_map_b: linear_table_b,
                gamma_map_r: gamma_table_r,
                gamma_map_g: gamma_table_g,
                gamma_map_b: gamma_table_b,
                im_stage: Some(im_stage),
                tone_map,
                params,
            }))
        }
        MappingColorSpace::Yrg(params) => {
            // We need to adapt PCS D50 to CIE XYZ 2006 with D65 white point first.
            let mut to_xyz = input_color_space.rgb_to_xyz_matrix().to_f32();
            to_xyz = to_xyz.mul_row::<1>(1.05785528f32);
            let output_d50 = output_color_space.rgb_to_xyz_matrix().to_f32();
            let mut to_rgb = output_d50.inverse();
            to_rgb = to_rgb.mul_row::<1>(1. / 1.05785528f32);

            Ok(Arc::new(ToneMapperImplYrg::<u16, 65536, CN, 65536> {
                linear_map_r: linear_table_r,
                linear_map_g: linear_table_g,
                linear_map_b: linear_table_b,
                gamma_map_r: gamma_table_r,
                gamma_map_g: gamma_table_g,
                gamma_map_b: gamma_table_b,
                to_xyz,
                to_rgb,
                tone_map,
                parameters: params,
            }))
        }
        MappingColorSpace::Jzazbz(brightness) => {
            // We need to adapt PCS D50 to XYZ with D65 white point first.
            let to_xyz = input_color_space.rgb_to_xyz_matrix().to_f32();
            let output_d65 = output_color_space.rgb_to_xyz_matrix().to_f32();
            let to_rgb = output_d65.inverse();

            Ok(Arc::new(ToneMapperImplJzazbz::<u16, 65536, CN, 65536> {
                linear_map_r: linear_table_r,
                linear_map_g: linear_table_g,
                linear_map_b: linear_table_b,
                gamma_map_r: gamma_table_r,
                gamma_map_g: gamma_table_g,
                gamma_map_b: gamma_table_b,
                to_xyz,
                to_rgb,
                tone_map,
                parameters: brightness,
            }))
        }
    }
}

#[cfg(all(target_arch = "aarch64", feature = "neon"))]
#[allow(clippy::borrowed_box)]
fn get_neon_transform<
    T: MangledCoercion + Clone + Copy + Default + Send + Sync + 'static + Debug,
    const N: usize,
    const CN: usize,
>(
    input_color_space: &ColorProfile,
    method: &ToneMappingMethod,
    working_color_space: &MappingColorSpace,
    linear_table_r: &Box<[f32; N]>,
    gamma_table_r: &Box<[T; 65536]>,
    tone_map: &Arc<SyncToneMap>,
    conversion: Matrix3f,
    bit_depth: usize,
    gamma_lut: usize,
) -> Option<Result<ToneMapping<T>, ForgeError>>
where
    u32: AsPrimitive<T>,
{
    if !matches!(working_color_space, MappingColorSpace::Rgb(_)) {
        return None;
    }
    let rgb_param = if let MappingColorSpace::Rgb(rgb_param) = &working_color_space
        && rgb_param.gamut_clipping == GamutClipping::NoClip
    {
        rgb_param
    } else {
        return None;
    };
    if let ToneMappingMethod::TunedReinhard(params) = &method {
        use crate::neon::{DisplayReinhardParamsNeon, HotTunedReinhardNeon};
        let primaries = input_color_space.rgb_to_xyz_matrix().to_f32();
        let luma_primaries: [f32; 3] = primaries.v[1];
        return Some(Ok(Arc::new(HotTunedReinhardNeon::<T, N, CN> {
            linear: linear_table_r.clone(),
            gamma: gamma_table_r.clone(),
            gamma_lut,
            adaptation_matrix: conversion,
            bit_depth,
            mapper: DisplayReinhardParamsNeon::new(
                params.content_max_brightness,
                params.display_max_brightness,
                203f32,
                luma_primaries,
            ),
            tone_map: tone_map.clone(),
            exposure: rgb_param.exposure,
        })));
    } else if let ToneMappingMethod::ExtendedReinhard { max_luma } = &method {
        use crate::neon::{ExtendedReinhardNeon, HotExtendedReinhardNeon};
        let primaries = input_color_space.rgb_to_xyz_matrix().to_f32();
        let luma_primaries: [f32; 3] = primaries.v[1];
        return Some(Ok(Arc::new(HotExtendedReinhardNeon::<T, N, CN> {
            linear: linear_table_r.clone(),
            gamma: gamma_table_r.clone(),
            gamma_lut,
            adaptation_matrix: conversion,
            bit_depth,
            mapper: ExtendedReinhardNeon::new(luma_primaries, *max_luma),
            tone_map: tone_map.clone(),
            exposure: rgb_param.exposure,
        })));
    } else if let ToneMappingMethod::ReinhardJodie = &method {
        use crate::neon::{HotReinhardJodieNeon, ReinhardJodieNeon};
        let primaries = input_color_space.rgb_to_xyz_matrix().to_f32();
        let luma_primaries: [f32; 3] = primaries.v[1];
        return Some(Ok(Arc::new(HotReinhardJodieNeon::<T, N, CN> {
            linear: linear_table_r.clone(),
            gamma: gamma_table_r.clone(),
            gamma_lut,
            adaptation_matrix: conversion,
            bit_depth,
            mapper: ReinhardJodieNeon::new(luma_primaries),
            tone_map: tone_map.clone(),
            exposure: rgb_param.exposure,
        })));
    } else if let ToneMappingMethod::Filmic = &method {
        use crate::neon::{HableNeon, HotFilmicNeon};
        return Some(Ok(Arc::new(HotFilmicNeon::<T, N, CN> {
            linear: linear_table_r.clone(),
            gamma: gamma_table_r.clone(),
            gamma_lut,
            adaptation_matrix: conversion,
            bit_depth,
            mapper: HableNeon,
            tone_map: tone_map.clone(),
            exposure: rgb_param.exposure,
        })));
    } else if let ToneMappingMethod::Aces = &method {
        use crate::neon::{AcesNeon, HotAcesNeon};
        return Some(Ok(Arc::new(HotAcesNeon::<T, N, CN> {
            linear: linear_table_r.clone(),
            gamma: gamma_table_r.clone(),
            gamma_lut,
            adaptation_matrix: conversion,
            bit_depth,
            mapper: AcesNeon::default(),
            tone_map: tone_map.clone(),
            exposure: rgb_param.exposure,
        })));
    } else if let ToneMappingMethod::Reinhard = &method {
        use crate::neon::{HotReinhardNeon, ReinhardNeon};
        return Some(Ok(Arc::new(HotReinhardNeon::<T, N, CN> {
            linear: linear_table_r.clone(),
            gamma: gamma_table_r.clone(),
            gamma_lut,
            adaptation_matrix: conversion,
            bit_depth,
            mapper: ReinhardNeon::default(),
            tone_map: tone_map.clone(),
            exposure: rgb_param.exposure,
        })));
    }
    None
}

#[cfg(all(target_arch = "x86_64", feature = "avx"))]
#[allow(clippy::borrowed_box)]
fn gen_avx_transform<
    T: MangledCoercion + Clone + Copy + Default + Send + Sync + 'static + Debug,
    const N: usize,
    const CN: usize,
>(
    input_color_space: &ColorProfile,
    method: &ToneMappingMethod,
    working_color_space: &MappingColorSpace,
    linear_table_r: &Box<[f32; N]>,
    gamma_table_r: &Box<[T; 65536]>,
    tone_map: &Arc<SyncToneMap>,
    conversion: Matrix3f,
    bit_depth: usize,
    gamma_lut: usize,
) -> Option<Result<ToneMapping<T>, ForgeError>>
where
    u32: AsPrimitive<T>,
{
    if !std::arch::is_x86_feature_detected!("avx2") || !std::arch::is_x86_feature_detected!("fma") {
        return None;
    }
    if !matches!(working_color_space, MappingColorSpace::Rgb(_)) {
        return None;
    }
    let rgb_param = if let MappingColorSpace::Rgb(rgb_param) = &working_color_space
        && rgb_param.gamut_clipping == GamutClipping::NoClip
    {
        rgb_param
    } else {
        return None;
    };
    if let ToneMappingMethod::ReinhardJodie = &method {
        use crate::avx::{HotReinhardJodieAvx, ReinhardJodieAvx};
        let primaries = input_color_space.rgb_to_xyz_matrix().to_f32();
        let luma_primaries: [f32; 3] = primaries.v[1];
        return Some(Ok(Arc::new(HotReinhardJodieAvx::<T, N, CN> {
            linear: linear_table_r.clone(),
            gamma: gamma_table_r.clone(),
            gamma_lut,
            adaptation_matrix: conversion,
            bit_depth,
            mapper: ReinhardJodieAvx::new(luma_primaries),
            tone_map: tone_map.clone(),
            exposure: rgb_param.exposure,
        })));
    } else if let ToneMappingMethod::Reinhard = &method {
        use crate::avx::{HotReinhardAvx, ReinhardAvx};
        return Some(Ok(Arc::new(HotReinhardAvx::<T, N, CN> {
            linear: linear_table_r.clone(),
            gamma: gamma_table_r.clone(),
            gamma_lut,
            adaptation_matrix: conversion,
            bit_depth,
            mapper: ReinhardAvx::default(),
            tone_map: tone_map.clone(),
            exposure: rgb_param.exposure,
        })));
    } else if let ToneMappingMethod::Filmic = &method {
        use crate::avx::{HableAvx, HotHableAvx};
        return Some(Ok(Arc::new(HotHableAvx::<T, N, CN> {
            linear: linear_table_r.clone(),
            gamma: gamma_table_r.clone(),
            gamma_lut,
            adaptation_matrix: conversion,
            bit_depth,
            mapper: HableAvx,
            tone_map: tone_map.clone(),
            exposure: rgb_param.exposure,
        })));
    } else if let ToneMappingMethod::Aces = &method {
        use crate::avx::{AcesAvx, HotAcesAvx};
        return Some(Ok(Arc::new(HotAcesAvx::<T, N, CN> {
            linear: linear_table_r.clone(),
            gamma: gamma_table_r.clone(),
            gamma_lut,
            adaptation_matrix: conversion,
            bit_depth,
            mapper: AcesAvx::default(),
            tone_map: tone_map.clone(),
            exposure: rgb_param.exposure,
        })));
    } else if let ToneMappingMethod::TunedReinhard(params) = &method {
        use crate::avx::{DisplayReinhardParamsAvx, HotTunedReinhardAvx};
        let primaries = input_color_space.rgb_to_xyz_matrix().to_f32();
        let luma_primaries: [f32; 3] = primaries.v[1];
        return Some(Ok(Arc::new(HotTunedReinhardAvx::<T, N, CN> {
            linear: linear_table_r.clone(),
            gamma: gamma_table_r.clone(),
            gamma_lut,
            adaptation_matrix: conversion,
            bit_depth,
            mapper: DisplayReinhardParamsAvx::new(
                params.content_max_brightness,
                params.display_max_brightness,
                203f32,
                luma_primaries,
            ),
            tone_map: tone_map.clone(),
            exposure: rgb_param.exposure,
        })));
    } else if let ToneMappingMethod::ExtendedReinhard { max_luma } = &method {
        use crate::avx::{ExtendedReinhardAvx, HotExtendedReinhardAvx};
        let primaries = input_color_space.rgb_to_xyz_matrix().to_f32();
        let luma_primaries: [f32; 3] = primaries.v[1];
        return Some(Ok(Arc::new(HotExtendedReinhardAvx::<T, N, CN> {
            linear: linear_table_r.clone(),
            gamma: gamma_table_r.clone(),
            gamma_lut,
            adaptation_matrix: conversion,
            bit_depth,
            mapper: ExtendedReinhardAvx::new(luma_primaries, *max_luma),
            tone_map: tone_map.clone(),
            exposure: rgb_param.exposure,
        })));
    }
    None
}

fn make_mapper<const CN: usize>(
    input_color_space: &ColorProfile,
    method: ToneMappingMethod,
) -> Arc<SyncToneMap> {
    let primaries = input_color_space.rgb_to_xyz_matrix().to_f32();
    let luma_primaries: [f32; 3] = primaries.v[1];
    let tone_map: Arc<SyncToneMap> = match method {
        ToneMappingMethod::Itu2408(data) => Arc::new(Rec2408ToneMapper::<CN>::new(
            data.content_max_brightness,
            data.display_max_brightness,
            luma_primaries,
        )),
        ToneMappingMethod::TunedReinhard(data) => {
            use crate::mappers::TunedReinhardToneMapper;
            Arc::new(TunedReinhardToneMapper::<CN>::new(
                data.content_max_brightness,
                data.display_max_brightness,
                203f32,
                luma_primaries,
            ))
        }
        ToneMappingMethod::Filmic => Arc::new(FilmicToneMapper::<CN>::default()),
        ToneMappingMethod::Aces => Arc::new(AcesToneMapper::<CN>::default()),
        ToneMappingMethod::ExtendedReinhard { max_luma } => {
            use crate::mappers::ExtendedReinhardToneMapper;
            Arc::new(ExtendedReinhardToneMapper::<CN> {
                primaries: luma_primaries,
                recip_max_l_sqr: 1. / (max_luma * max_luma),
            })
        }
        ToneMappingMethod::ReinhardJodie => Arc::new(ReinhardJodieToneMapper::<CN> {
            primaries: luma_primaries,
        }),
        ToneMappingMethod::Reinhard => Arc::new(ReinhardToneMapper::<CN>::default()),
        ToneMappingMethod::Clamp => Arc::new(ClampToneMapper::<CN>::default()),
        ToneMappingMethod::FilmicSpline(params) => {
            let spline = create_spline(params);
            Arc::new(SplineToneMapper::<CN> {
                spline,
                primaries: luma_primaries,
            })
        }
        ToneMappingMethod::Agx(look) => match look {
            AgxLook::Agx => Arc::new(AgxToneMapper::<CN> {
                primaries: luma_primaries,
                agx_custom_look: AgxDefault::custom_look(),
            }),
            AgxLook::Punchy => Arc::new(AgxToneMapper::<CN> {
                primaries: luma_primaries,
                agx_custom_look: AgxPunchy::custom_look(),
            }),
            AgxLook::Golden => Arc::new(AgxToneMapper::<CN> {
                primaries: luma_primaries,
                agx_custom_look: AgxGolden::custom_look(),
            }),
            AgxLook::Custom(look) => Arc::new(AgxToneMapper::<CN> {
                primaries: luma_primaries,
                agx_custom_look: look,
            }),
        },
    };
    tone_map
}

macro_rules! define8 {
    ($method: ident, $cn: expr, $name: expr) => {
        #[doc = concat!("Creates an ", $name," tone mapper. \
        \
        ICC profile do expect that for HDR tone management `CICP` tag will be used. \
        Tone mapper will search for `CICP` in [ColorProfile] and if there is some value, \
        then transfer function from `CICP` will be used. \
        Otherwise, we will interpolate ICC tone reproduction LUT tables.")]
        pub fn $method(
            input_color_space: &ColorProfile,
            output_color_space: &ColorProfile,
            method: ToneMappingMethod,
            working_color_space: MappingColorSpace,
        ) -> Result<ToneMapping<u8>, ForgeError> {
            create_tone_mapper_u8::<$cn>(
                input_color_space,
                output_color_space,
                method,
                working_color_space,
            )
        }
    };
}

define8!(create_tone_mapper_rgb, 3, "RGB8");
define8!(create_tone_mapper_rgba, 4, "RGBA8");

macro_rules! define16 {
    ($method: ident, $cn: expr, $bp: expr, $name: expr) => {
        #[doc = concat!("Creates an ", $name," tone mapper. \
        \
        ICC profile do expect that for HDR tone management `CICP` tag will be used. \
        Tone mapper will search for `CICP` in [ColorProfile] and if there is some value, \
        then transfer function from `CICP` will be used. \
        Otherwise, we will interpolate ICC tone reproduction LUT tables.")]
        pub fn $method(
            input_color_space: &ColorProfile,
            output_color_space: &ColorProfile,
            method: ToneMappingMethod,
            working_color_space: MappingColorSpace,
        ) -> Result<ToneMapping<u16>, ForgeError> {
            create_tone_mapper_u16::<$cn, $bp>(
                input_color_space,
                output_color_space,
                method,
                working_color_space,
            )
        }
    };
}

define16!(create_tone_mapper_rgb10, 3, 10, "RGB10");
define16!(create_tone_mapper_rgba10, 4, 10, "RGBA10");

define16!(create_tone_mapper_rgb12, 3, 12, "RGB12");
define16!(create_tone_mapper_rgba12, 4, 12, "RGBA12");

define16!(create_tone_mapper_rgb14, 3, 14, "RGB14");
define16!(create_tone_mapper_rgba14, 4, 14, "RGBA14");

define16!(create_tone_mapper_rgb16, 3, 16, "RGB16");
define16!(create_tone_mapper_rgba16, 4, 16, "RGBA16");