ojph_colour.cpp

Go to the documentation of this file.

//***************************************************************************/
// This software is released under the 2-Clause BSD license, included
// below.
//
// Copyright (c) 2019, Aous Naman
// Copyright (c) 2019, Kakadu Software Pty Ltd, Australia
// Copyright (c) 2019, The University of New South Wales, Australia
//
// 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.
//
// 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.
//***************************************************************************/
// This file is part of the OpenJPH software implementation.
// File: ojph_colour.cpp
// Author: Aous Naman
// Date: 28 August 2019
//***************************************************************************/
 
#include <cmath>
#include <climits>
#include <mutex>
 
#include "ojph_defs.h"
#include "ojph_arch.h"
#include "ojph_mem.h"
#include "ojph_colour.h"
#include "ojph_colour_local.h"
 
#include "ojph_params.h"
#include "../codestream/ojph_params_local.h"
 
namespace ojph {
 
  // defined elsewhere
  class line_buf;
 
  namespace local {

    void (*rev_convert)
      (const line_buf *src_line, const ui32 src_line_offset,
       line_buf *dst_line, const ui32 dst_line_offset,
       si64 shift, ui32 width) = NULL;

    void (*rev_encode_nlt)
      (const line_buf *src_line, const ui32 src_line_offset,
       line_buf *dst_line, ui32 bit_depth, bool is_signed, ui32 width,
       const nlt_rec *rec) = NULL;

    void (*rev_decode_nlt)
      (const line_buf *src_line, const ui32 src_line_offset,
       line_buf *dst_line, const ui32 dst_line_offset, ui32 bit_depth,
       bool is_signed, ui32 width, const nlt_rec *rec) = NULL;

    void (*rev_convert_nlt_type3)
      (const line_buf *src_line, const ui32 src_line_offset,
       line_buf *dst_line, const ui32 dst_line_offset,
       si64 shift, ui32 width) = NULL;

    void (*irv_convert_to_integer) (
      const line_buf *src_line, line_buf *dst_line, ui32 dst_line_offset,
      ui32 bit_depth, bool is_signed, ui32 width) = NULL;

    void (*irv_convert_to_float) (
      const line_buf *src_line, ui32 src_line_offset,
      line_buf *dst_line, ui32 bit_depth, bool is_signed, ui32 width) = NULL;

    void (*irv_convert_to_integer_nlt) (
      const line_buf *src_line, line_buf *dst_line, ui32 dst_line_offset,
      ui32 bit_depth, bool is_signed, ui32 width, const nlt_rec *rec) = NULL;

    void (*irv_convert_to_integer_nlt_type3) (
      const line_buf *src_line, line_buf *dst_line, ui32 dst_line_offset,
      ui32 bit_depth, bool is_signed, ui32 width) = NULL;

    void (*irv_convert_to_float_nlt) (
      const line_buf *src_line, ui32 src_line_offset, line_buf *dst_line,
      ui32 bit_depth, bool is_signed, ui32 width, const nlt_rec *rec) = NULL;

    void (*irv_convert_to_float_nlt_type3) (
      const line_buf *src_line, ui32 src_line_offset,
      line_buf *dst_line, ui32 bit_depth, bool is_signed, ui32 width) = NULL;

    void (*rct_forward)
      (const line_buf* r, const line_buf* g, const line_buf* b,
       line_buf* y, line_buf* cb, line_buf* cr, ui32 repeat) = NULL;

    void (*rct_backward)
      (const line_buf* r, const line_buf* g, const line_buf* b,
       line_buf* y, line_buf* cb, line_buf* cr, ui32 repeat) = NULL;

    void (*ict_forward)
      (const float *r, const float *g, const float *b,
       float *y, float *cb, float *cr, ui32 repeat) = NULL;

    void (*ict_backward)
      (const float *y, const float *cb, const float *cr,
       float *r, float *g, float *b, ui32 repeat) = NULL;

    void init_colour_transform_functions()
    {
      static std::once_flag colour_transform_functions_init_flag;
      std::call_once(colour_transform_functions_init_flag, []() {
#if !defined(OJPH_ENABLE_WASM_SIMD) || !defined(OJPH_EMSCRIPTEN)
 
        rev_convert = gen_rev_convert;
        rev_encode_nlt = gen_rev_encode_nlt;
        rev_decode_nlt = gen_rev_decode_nlt;
        rev_convert_nlt_type3 = gen_rev_convert_nlt_type3;
        irv_convert_to_integer = gen_irv_convert_to_integer;
        irv_convert_to_float = gen_irv_convert_to_float;
        irv_convert_to_integer_nlt = gen_irv_convert_to_integer_nlt;
        irv_convert_to_integer_nlt_type3 = gen_irv_convert_to_integer_nlt_type3;
        irv_convert_to_float_nlt = gen_irv_convert_to_float_nlt;
        irv_convert_to_float_nlt_type3 = gen_irv_convert_to_float_nlt_type3;
        rct_forward = gen_rct_forward;
        rct_backward = gen_rct_backward;
        ict_forward = gen_ict_forward;
        ict_backward = gen_ict_backward;
 
  #ifndef OJPH_DISABLE_SIMD
 
    #if (defined(OJPH_ARCH_X86_64) || defined(OJPH_ARCH_I386))
 
      #ifndef OJPH_DISABLE_SSE
        if (get_cpu_ext_level() >= X86_CPU_EXT_LEVEL_SSE)
        {
          ict_forward = sse_ict_forward;
          ict_backward = sse_ict_backward;
        }
      #endif // !OJPH_DISABLE_SSE
 
      #ifndef OJPH_DISABLE_SSE2
        if (get_cpu_ext_level() >= X86_CPU_EXT_LEVEL_SSE2)
        {
          rev_convert = sse2_rev_convert;
          rev_convert_nlt_type3 = sse2_rev_convert_nlt_type3;
          irv_convert_to_integer = sse2_irv_convert_to_integer;
          irv_convert_to_float = sse2_irv_convert_to_float;
          irv_convert_to_integer_nlt = sse2_irv_convert_to_integer_nlt;
          irv_convert_to_integer_nlt_type3 =
            sse2_irv_convert_to_integer_nlt_type3;
          irv_convert_to_float_nlt = sse2_irv_convert_to_float_nlt;
          irv_convert_to_float_nlt_type3 =
            sse2_irv_convert_to_float_nlt_type3;
          rct_forward = sse2_rct_forward;
          rct_backward = sse2_rct_backward;
        }
      #endif // !OJPH_DISABLE_SSE2
 
      #ifndef OJPH_DISABLE_AVX
        if (get_cpu_ext_level() >= X86_CPU_EXT_LEVEL_AVX)
        {
          ict_forward = avx_ict_forward;
          ict_backward = avx_ict_backward;
        }
      #endif // !OJPH_DISABLE_AVX
 
      #ifndef OJPH_DISABLE_AVX2
        if (get_cpu_ext_level() >= X86_CPU_EXT_LEVEL_AVX2)
        {
          rev_convert = avx2_rev_convert;
          rev_convert_nlt_type3 = avx2_rev_convert_nlt_type3;
          irv_convert_to_integer = avx2_irv_convert_to_integer;
          irv_convert_to_float = avx2_irv_convert_to_float;
          irv_convert_to_integer_nlt = avx2_irv_convert_to_integer_nlt;
          irv_convert_to_integer_nlt_type3 =
            avx2_irv_convert_to_integer_nlt_type3;
          irv_convert_to_float_nlt = avx2_irv_convert_to_float_nlt;
          irv_convert_to_float_nlt_type3 =
            avx2_irv_convert_to_float_nlt_type3;
          rct_forward = avx2_rct_forward;
          rct_backward = avx2_rct_backward;
        }
      #endif // !OJPH_DISABLE_AVX2
 
    #elif defined(OJPH_ARCH_ARM)
 
    #elif defined(OJPH_ARCH_PPC64LE)
 
        if (get_cpu_ext_level() >= PPC_CPU_EXT_LEVEL_ARCH_2_07)
        {
          // 128-bit VSX kernels; see ojph_simd_vsx.h
          rev_convert = vsx_rev_convert;
          rev_convert_nlt_type3 = vsx_rev_convert_nlt_type3;
          irv_convert_to_integer = vsx_irv_convert_to_integer;
          irv_convert_to_float = vsx_irv_convert_to_float;
          irv_convert_to_integer_nlt_type3 =
            vsx_irv_convert_to_integer_nlt_type3;
          irv_convert_to_float_nlt_type3 =
            vsx_irv_convert_to_float_nlt_type3;
          rct_forward = vsx_rct_forward;
          rct_backward = vsx_rct_backward;
          ict_forward = vsx_ict_forward;
          ict_backward = vsx_ict_backward;
        }
 
    #endif // !(defined(OJPH_ARCH_X86_64) || defined(OJPH_ARCH_I386))
 
  #endif // !OJPH_DISABLE_SIMD
 
#else // OJPH_ENABLE_WASM_SIMD
 
        rev_convert = wasm_rev_convert;
        rev_encode_nlt = gen_rev_encode_nlt;
        rev_decode_nlt = gen_rev_decode_nlt;
        rev_convert_nlt_type3 = wasm_rev_convert_nlt_type3;
        irv_convert_to_integer = wasm_irv_convert_to_integer;
        irv_convert_to_float = wasm_irv_convert_to_float;
        irv_convert_to_integer_nlt = gen_irv_convert_to_integer_nlt;
        irv_convert_to_integer_nlt_type3 = wasm_irv_convert_to_integer_nlt_type3;
        irv_convert_to_float_nlt = gen_irv_convert_to_float_nlt;
        irv_convert_to_float_nlt_type3 = wasm_irv_convert_to_float_nlt_type3;
        rct_forward = wasm_rct_forward;
        rct_backward = wasm_rct_backward;
        ict_forward = wasm_ict_forward;
        ict_backward = wasm_ict_backward;
 
#endif // !OJPH_ENABLE_WASM_SIMD
      });
    }

    const float CT_CNST::ALPHA_RF = 0.299f;
    const float CT_CNST::ALPHA_GF = 0.587f;
    const float CT_CNST::ALPHA_BF = 0.114f;
    const float CT_CNST::BETA_CbF = float(0.5/(1-double(CT_CNST::ALPHA_BF)));
    const float CT_CNST::BETA_CrF = float(0.5/(1-double(CT_CNST::ALPHA_RF)));
    const float CT_CNST::GAMMA_CB2G =
      float(2.0*double(ALPHA_BF)*(1.0-double(ALPHA_BF))/double(ALPHA_GF));
    const float CT_CNST::GAMMA_CR2G =
      float(2.0*double(ALPHA_RF)*(1.0-double(ALPHA_RF))/double(ALPHA_GF));
    const float CT_CNST::GAMMA_CB2B = float(2.0 * (1.0 - double(ALPHA_BF)));
    const float CT_CNST::GAMMA_CR2R = float(2.0 * (1.0 - double(ALPHA_RF)));


    template<int NLT_TYPE>
    static inline
    void local_gen_irv_convert_to_integer_nlt2or4(const line_buf *src_line,
      line_buf *dst_line, ui32 dst_line_offset,
      ui32 bit_depth, bool is_signed, ui32 width, const nlt_rec* rec)
    {
      assert((src_line->flags & line_buf::LFT_32BIT) &&
             (src_line->flags & line_buf::LFT_INTEGER) == 0 &&
             (dst_line->flags & line_buf::LFT_32BIT) &&
             (dst_line->flags & line_buf::LFT_INTEGER));
      ojph_unused(bit_depth);
      ojph_unused(is_signed);
 
      assert(rec->get_bit_depth() <= 32);
      const float* sp = src_line->f32;
      si32* dp = dst_line->i32 + dst_line_offset;
 
      float mul = rec->multiplier;
      float d_min = rec->fd_min;
      float d_max = rec->fd_max;
      float delta = rec->delta;
      float inv_delta = rec->inv_delta;
      float* lut = rec->dec_points;
 
      if (rec->is_signed())
      {
        const ui32 bias = (ui32)(1ULL << (rec->get_bit_depth() - 1)) + 1u;
        const float half = (float)(1ULL << (rec->get_bit_depth() - 1));
        for (int i = (int)width; i > 0; --i) {
          float t = *sp++ + 0.5f; // convert to [0, 1]
          t = ojph_max(t, d_min);
          t = ojph_min(t, d_max);
          ui32 k = (ui32)floorf((t - d_min) * inv_delta);
          float d_k = d_min + (float)k * delta;
          float t_k = lut[k];
          float t_kp1 = lut[k + 1];
          float z = t_k + (t - d_k) * inv_delta * (t_kp1 - t_k);
 
          si32 v = ojph_round(z * mul - half);
          if (NLT_TYPE == 4)
            v = (v >= 0) ? v : (si32)(-(ui32)v - bias);
          *dp++ = v;
        }
      }
      else
      {
        for (int i = (int)width; i > 0; --i) {
          float t = *sp++ + 0.5f;  // convert to [0, 1]
          t = ojph_max(t, d_min);
          t = ojph_min(t, d_max);
          ui32 k = (ui32)floorf((t - d_min) * inv_delta);
          float d_k = d_min + (float)k * delta;
          float t_k = lut[k];
          float t_kp1 = lut[k + 1];
          float z = t_k + (t - d_k) * inv_delta * (t_kp1 - t_k);
 
          si32 v = ojph_round(z * mul);
          *dp++ = v;
        }
      }
    }

    void gen_irv_convert_to_integer_nlt(const line_buf *src_line,
      line_buf *dst_line, ui32 dst_line_offset, ui32 bit_depth,
      bool is_signed, ui32 width, const nlt_rec* rec)
    {
      using nl = nlt_rec::nonlinearity;
      // if (rec->get_type() == nl::OJPH_NLT_NO_NLT)
      //   local_gen_irv_convert_to_integer_nlt0or3<false>(src_line, dst_line,
      //     dst_line_offset, bit_depth, is_signed, width);
      // else
      if (rec->get_type() == nl::OJPH_NLT_LUT_STYLE_NLT)
        local_gen_irv_convert_to_integer_nlt2or4<2>(src_line, dst_line,
          dst_line_offset, bit_depth, is_signed, width, rec);
      // else if (rec->get_type() == nl::OJPH_NLT_BINARY_COMPLEMENT_NLT)
      //   local_gen_irv_convert_to_integer_nlt0or3<true>(src_line, dst_line,
      //     dst_line_offset, bit_depth, is_signed, width);
      else if (rec->get_type() == nl::OJPH_NLT_BINARY_COMPLEMENT_PLUS_LUT)
        local_gen_irv_convert_to_integer_nlt2or4<4>(src_line, dst_line,
          dst_line_offset, bit_depth, is_signed, width, rec);
      else
        assert(0);
    }

    template<int NLT_TYPE>
    static inline
    void local_gen_irv_convert_to_float_nlt2or4(const line_buf *src_line,
      ui32 src_line_offset, line_buf *dst_line,
      ui32 bit_depth, bool is_signed, ui32 width, const nlt_rec* rec)
    {
      assert((src_line->flags & line_buf::LFT_32BIT) &&
             (src_line->flags & line_buf::LFT_INTEGER) &&
             (dst_line->flags & line_buf::LFT_32BIT) &&
             (dst_line->flags & line_buf::LFT_INTEGER) == 0);
      ojph_unused(is_signed);
      ojph_unused(bit_depth);
 
      assert(bit_depth <= 32);
      float mul = rec->multiplier;
      float ft_min = rec->ft_min;
      float ft_max = rec->ft_max;
      float delta = rec->delta;
      float inv_delta = rec->inv_delta;
      float* lut = rec->approx_enc_points;
 
      const si32* sp = src_line->i32 + src_line_offset;
      float* dp = dst_line->f32;
      if (rec->is_signed())
      {
        const ui32 bias = (ui32)(1ULL << (rec->get_bit_depth() - 1)) + 1u;
        for (int i = (int)width; i > 0; --i) {
          si32 v = *sp++;
          if (NLT_TYPE == 4)
            v = (v >= 0) ? v : (si32)(-(ui32)v - bias);
          float t = (float)v * mul + 0.5f;  // convert to [0, 1]
 
          t = ojph_max(t, ft_min);
          t = ojph_min(t, ft_max);
          ui32 k = (ui32)floorf((t - ft_min) * inv_delta);
          float d_k = ft_min + (float)k * delta;
          float t_k = lut[k];
          float t_kp1 = lut[k + 1];
          float y = t_k + (t - d_k) * inv_delta * (t_kp1 - t_k);
 
          *dp++ = y - 0.5f;
        }
      }
      else
      {
        for (int i = (int)width; i > 0; --i) {
          si32 v = *sp++;
          float t = (float)v * mul;  // it is in [0, 1]
 
          t = ojph_max(t, ft_min);
          t = ojph_min(t, ft_max);
          ui32 k = (ui32)floorf((t - ft_min) * inv_delta);
          float d_k = ft_min + (float)k * delta;
          float t_k = lut[k];
          float t_kp1 = lut[k + 1];
          float y = t_k + (t - d_k) * inv_delta * (t_kp1 - t_k);
 
          *dp++ = y - 0.5f;
        }
      }
    }

    // The precise inverse of the look-up table, for a table whose derived
    // version is fine enough that one step of it crosses at most one entry of
    // the table that was supplied.  The entry that a value falls in is found
    // with the same arithmetic that the approximation above uses, and the
    // entry after that is looked at as well, to catch the values that fall
    // past the entry that was found.
    template<int NLT_TYPE>
    static inline
    void local_gen_irv_convert_to_float_nlt2or4_precise_1step(
      const line_buf *src_line, ui32 src_line_offset, line_buf *dst_line,
      ui32 bit_depth, bool is_signed, ui32 width, const nlt_rec* rec)
    {
      assert((src_line->flags & line_buf::LFT_32BIT) &&
             (src_line->flags & line_buf::LFT_INTEGER) &&
             (dst_line->flags & line_buf::LFT_32BIT) &&
             (dst_line->flags & line_buf::LFT_INTEGER) == 0);
      ojph_unused(is_signed);
      ojph_unused(bit_depth);
 
      assert(bit_depth <= 32);
      float mul = rec->multiplier;
      float ft_min = rec->ft_min;
      float ft_max = rec->ft_max;
      float fd_min = rec->fd_min;
      float delta = rec->delta;
      float inv_delta = rec->inv_delta;
      ui16* indices= rec->precise_enc_dec_indices;
      float* lut = rec->precise_dec_lut;
 
      const si32* sp = src_line->i32 + src_line_offset;
      float* dp = dst_line->f32;
      if (rec->is_signed())
      {
        const ui32 bias = (ui32)(1ULL << (rec->get_bit_depth() - 1)) + 1u;
        for (int i = (int)width; i > 0; --i) {
          si32 v = *sp++;
          if (NLT_TYPE == 4)
            v = (v >= 0) ? v : (si32)(-(ui32)v - bias);
          float t = (float)v * mul + 0.5f;  // convert to [0, 1]
 
          t = ojph_max(t, ft_min);
          t = ojph_min(t, ft_max);
          ui32 k = (ui32)floorf((t - ft_min) * inv_delta);
          k = indices[k];
          float y0 = lut[k];      // lookup 0
          float y1 = lut[k + 1];  // lookup 1
          float y2 = lut[k + 2];  // lookup 2: one-step search
          if (t > y1)
          { y0 = y1; y1 = y2; ++k; }
          float y = fd_min + (float)k * delta +
            (y1 > y0 ? (t - y0) * delta / (y1 - y0) : 0.0f);
 
          *dp++ = y - 0.5f;
        }
      }
      else
      {
        for (int i = (int)width; i > 0; --i) {
          si32 v = *sp++;
          float t = (float)v * mul;  // it is in [0, 1]
 
          t = ojph_max(t, ft_min);
          t = ojph_min(t, ft_max);
          ui32 k = (ui32)floorf((t - ft_min) * inv_delta);
          k = indices[k];
          float y0 = lut[k];      // lookup 0
          float y1 = lut[k + 1];  // lookup 1
          float y2 = lut[k + 2];  // lookup 2: one-step search
          if (t > y1)
          { y0 = y1; y1 = y2; ++k; }
          float y = fd_min + (float)k * delta +
            (y1 > y0 ? (t - y0) * delta / (y1 - y0) : 0.0f);
 
          *dp++ = y - 0.5f;
        }
      }
    }

    // The precise inverse of the look-up table, for a table whose derived
    // version is coarse enough that one step of it can cross more than one
    // entry of the table that was supplied.  The entry that a value falls in
    // then has to be found by searching, which does not vectorize.  This is
    // the implementation that the SIMD versions in ojph_colour_sse2.cpp and
    // ojph_colour_avx2.cpp copy, so keep those in step when this changes.
    template<int NLT_TYPE>
    static inline
    void local_gen_irv_convert_to_float_nlt2or4_precise_search(
      const line_buf *src_line, ui32 src_line_offset, line_buf *dst_line,
      ui32 bit_depth, bool is_signed, ui32 width, const nlt_rec* rec)
    {
      assert((src_line->flags & line_buf::LFT_32BIT) &&
             (src_line->flags & line_buf::LFT_INTEGER) &&
             (dst_line->flags & line_buf::LFT_32BIT) &&
             (dst_line->flags & line_buf::LFT_INTEGER) == 0);
      ojph_unused(is_signed);
      ojph_unused(bit_depth);
 
      assert(bit_depth <= 32);
      assert(rec->precise_max_steps > 1);
      float mul = rec->multiplier;
      float ft_min = rec->ft_min;
      float ft_max = rec->ft_max;
      float fd_min = rec->fd_min;
      float delta = rec->delta;
      float inv_delta = rec->inv_delta;
      ui16* indices= rec->precise_enc_dec_indices;
      float* lut = rec->precise_dec_lut;
      ui32 num_points = rec->num_points;
 
      const si32* sp = src_line->i32 + src_line_offset;
      float* dp = dst_line->f32;
      if (rec->is_signed())
      {
        const ui32 bias = (ui32)(1ULL << (rec->get_bit_depth() - 1)) + 1u;
        for (int i = (int)width; i > 0; --i) {
          si32 v = *sp++;
          if (NLT_TYPE == 4)
            v = (v >= 0) ? v : (si32)(-(ui32)v - bias);
          float t = (float)v * mul + 0.5f;  // convert to [0, 1]
 
          t = ojph_max(t, ft_min);
          t = ojph_min(t, ft_max);
          ui32 k = (ui32)floorf((t - ft_min) * inv_delta);
          k = indices[k];
 
          ui32 lo = k;
          ui32 hi = ojph_min(k + rec->precise_max_steps + 1, num_points);
          while (hi - lo > 1)
          {
            ui32 mid = (lo + hi) >> 1;
            if (t >= lut[mid])
              lo = mid;
            else
              hi = mid;
          }
 
          float y0 = lut[lo];
          float y1 = lut[hi];
          float y = fd_min + (float)lo * delta +
            (y1 > y0 ? (t - y0) * delta / (y1 - y0) : 0.0f);
 
          *dp++ = y - 0.5f;
        }
      }
      else
      {
        for (int i = (int)width; i > 0; --i) {
          si32 v = *sp++;
          float t = (float)v * mul;  // it is in [0, 1]
 
          t = ojph_max(t, ft_min);
          t = ojph_min(t, ft_max);
          ui32 k = (ui32)floorf((t - ft_min) * inv_delta);
          k = indices[k];
 
          ui32 lo = k;
          ui32 hi = ojph_min(k + rec->precise_max_steps + 1, num_points);
          while (hi - lo > 1)
          {
            ui32 mid = (lo + hi) >> 1;
            if (t >= lut[mid])
              lo = mid;
            else
              hi = mid;
          }
 
          float y0 = lut[lo];
          float y1 = lut[hi];
          float y = fd_min + (float)lo * delta +
            (y1 > y0 ? (t - y0) * delta / (y1 - y0) : 0.0f);
 
          *dp++ = y - 0.5f;
        }
      }
    }

    void gen_irv_convert_to_float_nlt(const line_buf *src_line,
      ui32 src_line_offset, line_buf *dst_line, ui32 bit_depth,
      bool is_signed, ui32 width, const nlt_rec* rec)
    {
      using nl = nlt_rec::nonlinearity;
      // if (rec->get_type() == nl::OJPH_NLT_NO_NLT)
      //   local_gen_irv_convert_to_float_nlt0or3<false>(src_line,
      //     src_line_offset, dst_line, bit_depth, is_signed, width);
      // else
      if (rec->get_type() == nl::OJPH_NLT_LUT_STYLE_NLT)
      {
        if (rec->approx_encoding_nlt)
          local_gen_irv_convert_to_float_nlt2or4<2>(src_line,
            src_line_offset, dst_line, bit_depth, is_signed, width, rec);
        else if (rec->precise_max_steps == 1)
          local_gen_irv_convert_to_float_nlt2or4_precise_1step<2>(src_line,
            src_line_offset, dst_line, bit_depth, is_signed, width, rec);
        else
          local_gen_irv_convert_to_float_nlt2or4_precise_search<2>(src_line,
            src_line_offset, dst_line, bit_depth, is_signed, width, rec);
      }
      // else if (rec->get_type() == nl::OJPH_NLT_BINARY_COMPLEMENT_NLT)
      //   local_gen_irv_convert_to_float_nlt0or3<true>(src_line,
      //     src_line_offset, dst_line, bit_depth, is_signed, width);
      else if (rec->get_type() == nl::OJPH_NLT_BINARY_COMPLEMENT_PLUS_LUT)
      {
        if (rec->approx_encoding_nlt)
          local_gen_irv_convert_to_float_nlt2or4<4>(src_line,
            src_line_offset, dst_line, bit_depth, is_signed, width, rec);
        else if (rec->precise_max_steps == 1)
          local_gen_irv_convert_to_float_nlt2or4_precise_1step<4>(src_line,
            src_line_offset, dst_line, bit_depth, is_signed, width, rec);
        else
          local_gen_irv_convert_to_float_nlt2or4_precise_search<4>(src_line,
            src_line_offset, dst_line, bit_depth, is_signed, width, rec);
      }
      else
        assert(0);
    }

    // The LUT style nonlinearities, types 2 and 4, are implemented for the
    // irreversible wavelet only.  param_nlt::check_validity refuses to write a
    // codestream that combines them with the reversible wavelet, so this entry
    // point is unreachable while that check is in place.
    void gen_rev_encode_nlt(
      const line_buf *src_line, const ui32 src_line_offset,
      line_buf *dst_line, ui32 bit_depth, bool is_signed, ui32 width,
      const nlt_rec* rec)
    {
      ojph_unused(src_line);
      ojph_unused(src_line_offset);
      ojph_unused(dst_line);
      ojph_unused(bit_depth);
      ojph_unused(is_signed);
      ojph_unused(width);
      ojph_unused(rec);
      assert(0);
    }

    void gen_rev_decode_nlt(
      const line_buf *src_line, const ui32 src_line_offset,
      line_buf *dst_line, const ui32 dst_line_offset,
      ui32 bit_depth, bool is_signed, ui32 width, const nlt_rec* rec)
    {
      // A codestream can still carry a LUT style nonlinearity together with the
      // reversible wavelet, because other encoders may produce one; this
      // library cannot invert such a nonlinearity, so the samples are passed
      // through unchanged.  The caller was warned about that when the headers
      // of the codestream were read, which is why this is not a failure here.
      ojph_unused(src_line);
      ojph_unused(src_line_offset);
      ojph_unused(dst_line);
      ojph_unused(dst_line_offset);
      ojph_unused(bit_depth);
      ojph_unused(is_signed);
      ojph_unused(width);
      ojph_unused(rec);
      assert(0);
    }
 
#if !defined(OJPH_ENABLE_WASM_SIMD) || !defined(OJPH_EMSCRIPTEN)

    void gen_rev_convert(
      const line_buf *src_line, const ui32 src_line_offset,
      line_buf *dst_line, const ui32 dst_line_offset,
      si64 shift, ui32 width)
    {
      if (src_line->flags & line_buf::LFT_32BIT)
      {
        if (dst_line->flags & line_buf::LFT_32BIT)
        {
          const si32 *sp = src_line->i32 + src_line_offset;
          si32 *dp = dst_line->i32 + dst_line_offset;
          si32 s = (si32)shift;
          for (ui32 i = width; i > 0; --i)
            *dp++ = *sp++ + s;
        }
        else
        {
          const si32 *sp = src_line->i32 + src_line_offset;
          si64 *dp = dst_line->i64 + dst_line_offset;
          for (ui32 i = width; i > 0; --i)
            *dp++ = *sp++ + shift;
        }
      }
      else
      {
        assert(src_line->flags & line_buf::LFT_64BIT);
        assert(dst_line->flags & line_buf::LFT_32BIT);
        const si64 *sp = src_line->i64 + src_line_offset;
        si32 *dp = dst_line->i32 + dst_line_offset;
        for (ui32 i = width; i > 0; --i)
          *dp++ = (si32)(*sp++ + shift);
      }
    }

    void gen_rev_convert_nlt_type3(
      const line_buf *src_line, const ui32 src_line_offset,
      line_buf *dst_line, const ui32 dst_line_offset,
      si64 shift, ui32 width)
    {
      if (src_line->flags & line_buf::LFT_32BIT)
      {
        if (dst_line->flags & line_buf::LFT_32BIT)
        {
          const si32 *sp = src_line->i32 + src_line_offset;
          si32 *dp = dst_line->i32 + dst_line_offset;
          si32 s = (si32)shift;
          for (ui32 i = width; i > 0; --i) {
            const si32 v = *sp++;
            *dp++ = v >= 0 ? v : (- v - s);
          }
        }
        else
        {
          const si32 *sp = src_line->i32 + src_line_offset;
          si64 *dp = dst_line->i64 + dst_line_offset;
          for (ui32 i = width; i > 0; --i) {
            const si64 v = *sp++;
            *dp++ = v >= 0 ? v : (- v - shift);
          }
        }
      }
      else
      {
        assert(src_line->flags & line_buf::LFT_64BIT);
        assert(dst_line->flags & line_buf::LFT_32BIT);
        const si64 *sp = src_line->i64 + src_line_offset;
        si32 *dp = dst_line->i32 + dst_line_offset;
        for (ui32 i = width; i > 0; --i) {
          const si64 v = *sp++;
          *dp++ = (si32)(v >= 0 ? v : (- v - shift));
        }
      }
    }

    template<bool NLT_TYPE3>
    static inline
    void local_gen_irv_convert_to_integer_nlt0or3(const line_buf *src_line,
      line_buf *dst_line, ui32 dst_line_offset,
      ui32 bit_depth, bool is_signed, ui32 width)
    {
      assert((src_line->flags & line_buf::LFT_32BIT) &&
             (src_line->flags & line_buf::LFT_INTEGER) == 0 &&
             (dst_line->flags & line_buf::LFT_32BIT) &&
             (dst_line->flags & line_buf::LFT_INTEGER));
 
      assert(bit_depth <= 32);
      const float* sp = src_line->f32;
      si32* dp = dst_line->i32 + dst_line_offset;
      // There is the possibility that converting to integer will
      // exceed the dynamic range of 32bit integer; therefore, care must be
      // exercised.
      // We look if the floating point number is outside the half-closed
      // interval [-0.5f, 0.5f). If so, we limit the resulting integer
      // to the maximum/minimum that number supports.
      si32 neg_limit = (si32)INT_MIN >> (32 - bit_depth);
      float mul = (float)(1ull << bit_depth);
      float fl_up_lim = -(float)neg_limit; // val < upper
      float fl_low_lim = (float)neg_limit; // val >= lower
      si32 s32_up_lim = INT_MAX >> (32 - bit_depth);
      si32 s32_low_lim = INT_MIN >> (32 - bit_depth);
 
      if (is_signed)
      {
        const ui32 bias = (ui32)(1ULL << (bit_depth - 1)) + 1u;
        for (int i = (int)width; i > 0; --i) {
          float t = *sp++ * mul;
          si32 v = ojph_round(t);
          v = t >= fl_low_lim ? v : s32_low_lim;
          v = t <  fl_up_lim  ? v : s32_up_lim;
          if (NLT_TYPE3)
            v = (v >= 0) ? v : (si32)(-(ui32)v - bias);
          *dp++ = v;
        }
      }
      else
      {
        const si32 half = (si32)(1ULL << (bit_depth - 1));
        for (int i = (int)width; i > 0; --i) {
          float t = *sp++ * mul;
          si32 v = ojph_round(t);
          v = t >= fl_low_lim ? v : s32_low_lim;
          v = t <  fl_up_lim  ? v : s32_up_lim;
          *dp++ = v + half;
        }
      }
    }

    void gen_irv_convert_to_integer(const line_buf *src_line,
      line_buf *dst_line, ui32 dst_line_offset,
      ui32 bit_depth, bool is_signed, ui32 width)
    {
      local_gen_irv_convert_to_integer_nlt0or3<false>(src_line, dst_line,
        dst_line_offset, bit_depth, is_signed, width);
    }

    void gen_irv_convert_to_integer_nlt_type3(const line_buf *src_line,
      line_buf *dst_line, ui32 dst_line_offset,
      ui32 bit_depth, bool is_signed, ui32 width)
    {
      local_gen_irv_convert_to_integer_nlt0or3<true>(src_line, dst_line,
        dst_line_offset, bit_depth, is_signed, width);
    }

    template<bool NLT_TYPE3>
    static inline
    void local_gen_irv_convert_to_float_nlt0or3(const line_buf *src_line,
      ui32 src_line_offset, line_buf *dst_line,
      ui32 bit_depth, bool is_signed, ui32 width)
    {
      assert((src_line->flags & line_buf::LFT_32BIT) &&
             (src_line->flags & line_buf::LFT_INTEGER) &&
             (dst_line->flags & line_buf::LFT_32BIT) &&
             (dst_line->flags & line_buf::LFT_INTEGER) == 0);
 
      assert(bit_depth <= 32);
      float mul = (float)(1.0 / (double)(1ULL << bit_depth));
 
      const si32* sp = src_line->i32 + src_line_offset;
      float* dp = dst_line->f32;
      if (is_signed)
      {
        const ui32 bias = (ui32)(1ULL << (bit_depth - 1)) + 1u;
        for (int i = (int)width; i > 0; --i) {
          si32 v = *sp++;
          if (NLT_TYPE3)
            v = (v >= 0) ? v : (si32)(-(ui32)v - bias);
          *dp++ = (float)v * mul;
        }
      }
      else
      {
        const si32 half = (si32)(1ULL << (bit_depth - 1));
        for (int i = (int)width; i > 0; --i) {
          si32 v = *sp++;
          v -= half;
          *dp++ = (float)v * mul;
        }
      }
    }

    void gen_irv_convert_to_float(const line_buf *src_line,
      ui32 src_line_offset, line_buf *dst_line,
      ui32 bit_depth, bool is_signed, ui32 width)
    {
      local_gen_irv_convert_to_float_nlt0or3<false>(src_line, src_line_offset,
        dst_line, bit_depth, is_signed, width);
    }

    void gen_irv_convert_to_float_nlt_type3(const line_buf *src_line,
      ui32 src_line_offset, line_buf *dst_line,
      ui32 bit_depth, bool is_signed, ui32 width)
    {
      local_gen_irv_convert_to_float_nlt0or3<true>(src_line, src_line_offset,
        dst_line, bit_depth, is_signed, width);
    }

    void gen_rct_forward(
      const line_buf *r, const line_buf *g, const line_buf *b,
      line_buf *y, line_buf *cb, line_buf *cr, ui32 repeat)
    {
      assert((y->flags  & line_buf::LFT_INTEGER) &&
             (cb->flags & line_buf::LFT_INTEGER) &&
             (cr->flags & line_buf::LFT_INTEGER) &&
             (r->flags  & line_buf::LFT_INTEGER) &&
             (g->flags  & line_buf::LFT_INTEGER) &&
             (b->flags  & line_buf::LFT_INTEGER));
 
      if  (y->flags & line_buf::LFT_32BIT)
      {
        assert((y->flags  & line_buf::LFT_32BIT) &&
               (cb->flags & line_buf::LFT_32BIT) &&
               (cr->flags & line_buf::LFT_32BIT) &&
               (r->flags  & line_buf::LFT_32BIT) &&
               (g->flags  & line_buf::LFT_32BIT) &&
               (b->flags  & line_buf::LFT_32BIT));
        const si32 *rp = r->i32, * gp = g->i32, * bp = b->i32;
        si32 *yp = y->i32, * cbp = cb->i32, * crp = cr->i32;
        for (ui32 i = repeat; i > 0; --i)
        {
          si32 rr = *rp++, gg = *gp++, bb = *bp++;
          *yp++ = (rr + 2 * gg + bb) >> 2;
          *cbp++ = (bb - gg);
          *crp++ = (rr - gg);
        }
      }
      else
      {
        assert((y->flags  & line_buf::LFT_64BIT) &&
               (cb->flags & line_buf::LFT_64BIT) &&
               (cr->flags & line_buf::LFT_64BIT) &&
               (r->flags  & line_buf::LFT_32BIT) &&
               (g->flags  & line_buf::LFT_32BIT) &&
               (b->flags  & line_buf::LFT_32BIT));
        const si32 *rp = r->i32, *gp = g->i32, *bp = b->i32;
        si64 *yp = y->i64, *cbp = cb->i64, *crp = cr->i64;
        for (ui32 i = repeat; i > 0; --i)
        {
          si64 rr = *rp++, gg = *gp++, bb = *bp++;
          *yp++ = (rr + 2 * gg + bb) >> 2;
          *cbp++ = (bb - gg);
          *crp++ = (rr - gg);
        }
      }
    }

    void gen_rct_backward(
      const line_buf *y, const line_buf *cb, const line_buf *cr,
      line_buf *r, line_buf *g, line_buf *b, ui32 repeat)
    {
      assert((y->flags  & line_buf::LFT_INTEGER) &&
             (cb->flags & line_buf::LFT_INTEGER) &&
             (cr->flags & line_buf::LFT_INTEGER) &&
             (r->flags  & line_buf::LFT_INTEGER) &&
             (g->flags  & line_buf::LFT_INTEGER) &&
             (b->flags  & line_buf::LFT_INTEGER));
 
      if (y->flags & line_buf::LFT_32BIT)
      {
        assert((y->flags  & line_buf::LFT_32BIT) &&
               (cb->flags & line_buf::LFT_32BIT) &&
               (cr->flags & line_buf::LFT_32BIT) &&
               (r->flags  & line_buf::LFT_32BIT) &&
               (g->flags  & line_buf::LFT_32BIT) &&
               (b->flags  & line_buf::LFT_32BIT));
        const si32 *yp = y->i32, *cbp = cb->i32, *crp = cr->i32;
        si32 *rp = r->i32, *gp = g->i32, *bp = b->i32;
        for (ui32 i = repeat; i > 0; --i)
        {
          si32 yy = *yp++, cbb = *cbp++, crr = *crp++;
          si32 gg = yy - ((cbb + crr) >> 2);
          *rp++ = crr + gg;
          *gp++ = gg;
          *bp++ = cbb + gg;
        }
      }
      else
      {
        assert((y->flags  & line_buf::LFT_64BIT) &&
               (cb->flags & line_buf::LFT_64BIT) &&
               (cr->flags & line_buf::LFT_64BIT) &&
               (r->flags  & line_buf::LFT_32BIT) &&
               (g->flags  & line_buf::LFT_32BIT) &&
               (b->flags  & line_buf::LFT_32BIT));
        const si64 *yp = y->i64, *cbp = cb->i64, *crp = cr->i64;
        si32 *rp = r->i32, *gp = g->i32, *bp = b->i32;
        for (ui32 i = repeat; i > 0; --i)
        {
          si64 yy = *yp++, cbb = *cbp++, crr = *crp++;
          si64 gg = yy - ((cbb + crr) >> 2);
          *rp++ = (si32)(crr + gg);
          *gp++ = (si32)gg;
          *bp++ = (si32)(cbb + gg);
        }
      }
    }

    void gen_ict_forward(const float *r, const float *g, const float *b,
                         float *y, float *cb, float *cr, ui32 repeat)
    {
      for (ui32 i = repeat; i > 0; --i)
      {
        *y = CT_CNST::ALPHA_RF * *r
           + CT_CNST::ALPHA_GF * *g++
           + CT_CNST::ALPHA_BF * *b;
        *cb++ = CT_CNST::BETA_CbF * (*b++ - *y);
        *cr++ = CT_CNST::BETA_CrF * (*r++ - *y++);
      }
    }

    void gen_ict_backward(const float *y, const float *cb, const float *cr,
                          float *r, float *g, float *b, ui32 repeat)
    {
      for (ui32 i = repeat; i > 0; --i)
      {
        *g++ = *y - CT_CNST::GAMMA_CR2G * *cr - CT_CNST::GAMMA_CB2G * *cb;
        *r++ = *y + CT_CNST::GAMMA_CR2R * *cr++;
        *b++ = *y++ + CT_CNST::GAMMA_CB2B * *cb++;
      }
    }
 
#endif // !OJPH_ENABLE_WASM_SIMD
 
  }
}