ojph_transform.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_transform.cpp
// Author: Aous Naman
// Date: 28 August 2019
//***************************************************************************/
 
#include <cstdio>
#include <mutex>
 
#include "ojph_arch.h"
#include "ojph_mem.h"
#include "ojph_transform.h"
#include "ojph_transform_local.h"
#include "ojph_params.h"
#include "../codestream/ojph_params_local.h"
 
namespace ojph {
 
  // defined elsewhere
  class line_buf;
 
  namespace local {

    // Reversible functions

    void (*rev_vert_step)
      (const lifting_step* s, const line_buf* sig, const line_buf* other,
        const line_buf* aug, ui32 repeat, bool synthesis) = NULL;

    void (*rev_horz_ana)
      (const param_atk* atk, const line_buf* ldst, const line_buf* hdst,
        const line_buf* src, ui32 width, bool even) = NULL;

    void (*rev_horz_syn)
      (const param_atk* atk, const line_buf* dst, const line_buf* lsrc,
        const line_buf* hsrc, ui32 width, bool even) = NULL;
    
    // Irreversible functions

    void (*irv_vert_step)
      (const lifting_step* s, const line_buf* sig, const line_buf* other,
        const line_buf* aug, ui32 repeat, bool synthesis) = NULL;

    void (*irv_vert_times_K)
      (float K, const line_buf* aug, ui32 repeat) = NULL;

    void (*irv_horz_ana)
      (const param_atk* atk, const line_buf* ldst, const line_buf* hdst,
        const line_buf* src, ui32 width, bool even) = NULL;

    void (*irv_horz_syn)
      (const param_atk* atk, const line_buf* dst, const line_buf* lsrc,
        const line_buf* hsrc, ui32 width, bool even) = NULL;

    void init_wavelet_transform_functions()
    {
      static std::once_flag wavelet_transform_functions_init_flag;
      std::call_once(wavelet_transform_functions_init_flag, [](){
#if !defined(OJPH_ENABLE_WASM_SIMD) || !defined(OJPH_EMSCRIPTEN)
 
        rev_vert_step             = gen_rev_vert_step;
        rev_horz_ana              = gen_rev_horz_ana;
        rev_horz_syn              = gen_rev_horz_syn;
 
        irv_vert_step             = gen_irv_vert_step;
        irv_vert_times_K          = gen_irv_vert_times_K;
        irv_horz_ana              = gen_irv_horz_ana;
        irv_horz_syn              = gen_irv_horz_syn;
 
  #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)
        {
          irv_vert_step             = sse_irv_vert_step;
          irv_vert_times_K          = sse_irv_vert_times_K;
          irv_horz_ana              = sse_irv_horz_ana;
          irv_horz_syn              = sse_irv_horz_syn;
        }
      #endif // !OJPH_DISABLE_SSE
 
      #ifndef OJPH_DISABLE_SSE2
        if (get_cpu_ext_level() >= X86_CPU_EXT_LEVEL_SSE2)
        {
          rev_vert_step             = sse2_rev_vert_step;
          rev_horz_ana              = sse2_rev_horz_ana;
          rev_horz_syn              = sse2_rev_horz_syn;
        }
      #endif // !OJPH_DISABLE_SSE2
 
      #ifndef OJPH_DISABLE_AVX
        if (get_cpu_ext_level() >= X86_CPU_EXT_LEVEL_AVX)
        {
          irv_vert_step             = avx_irv_vert_step;
          irv_vert_times_K          = avx_irv_vert_times_K;
          irv_horz_ana              = avx_irv_horz_ana;      
          irv_horz_syn              = avx_irv_horz_syn;
        }
      #endif // !OJPH_DISABLE_AVX
 
      #ifndef OJPH_DISABLE_AVX2
        if (get_cpu_ext_level() >= X86_CPU_EXT_LEVEL_AVX2)
        {
          rev_vert_step             = avx2_rev_vert_step;
          rev_horz_ana              = avx2_rev_horz_ana;
          rev_horz_syn              = avx2_rev_horz_syn;
        }
      #endif // !OJPH_DISABLE_AVX2
 
      #if (defined(OJPH_ARCH_X86_64) && !defined(OJPH_DISABLE_AVX512))
        if (get_cpu_ext_level() >= X86_CPU_EXT_LEVEL_AVX512)
        {
          // rev_vert_step             = avx512_rev_vert_step;
          // rev_horz_ana              = avx512_rev_horz_ana;
          // rev_horz_syn              = avx512_rev_horz_syn;
 
          irv_vert_step             = avx512_irv_vert_step;
          irv_vert_times_K          = avx512_irv_vert_times_K;
          irv_horz_ana              = avx512_irv_horz_ana;
          irv_horz_syn              = avx512_irv_horz_syn;
        }
      #endif // !OJPH_DISABLE_AVX512
    
    #elif defined(OJPH_ARCH_ARM)
 
    #endif // !(defined(OJPH_ARCH_X86_64) || defined(OJPH_ARCH_I386))
 
  #endif // !OJPH_DISABLE_SIMD
 
#else // OJPH_ENABLE_WASM_SIMD
        rev_vert_step             = wasm_rev_vert_step;
        rev_horz_ana              = wasm_rev_horz_ana;
        rev_horz_syn              = wasm_rev_horz_syn;
        
        irv_vert_step             = wasm_irv_vert_step;
        irv_vert_times_K          = wasm_irv_vert_times_K;
        irv_horz_ana              = wasm_irv_horz_ana;
        irv_horz_syn              = wasm_irv_horz_syn;
#endif // !OJPH_ENABLE_WASM_SIMD
      });
    }
    
 
#if !defined(OJPH_ENABLE_WASM_SIMD) || !defined(OJPH_EMSCRIPTEN)

    static
    void gen_rev_vert_step32(const lifting_step* s, const line_buf* sig, 
                             const line_buf* other, const line_buf* aug, 
                             ui32 repeat, bool synthesis)
    {
      const si32 a = s->rev.Aatk;
      const si32 b = s->rev.Batk;
      const ui8 e = s->rev.Eatk;
 
      si32* dst = aug->i32;
      const si32* src1 = sig->i32, * src2 = other->i32;
      // The general definition of the wavelet in Part 2 is slightly 
      // different to part 2, although they are mathematically equivalent
      // here, we identify the simpler form from Part 1 and employ them
      if (a == 1)
      { // 5/3 update and any case with a == 1
        if (synthesis)
          for (ui32 i = repeat; i > 0; --i)
            *dst++ -= (b + *src1++ + *src2++) >> e;
        else
          for (ui32 i = repeat; i > 0; --i)
            *dst++ += (b + *src1++ + *src2++) >> e;
      }
      else if (a == -1 && b == 1 && e == 1)
      { // 5/3 predict
        if (synthesis)
          for (ui32 i = repeat; i > 0; --i)
            *dst++ += (*src1++ + *src2++) >> e;
        else
          for (ui32 i = repeat; i > 0; --i)
            *dst++ -= (*src1++ + *src2++) >> e;
      }
      else if (a == -1)
      { // any case with a == -1, which is not 5/3 predict
        if (synthesis)
          for (ui32 i = repeat; i > 0; --i)
            *dst++ -= (b - (*src1++ + *src2++)) >> e;
        else
          for (ui32 i = repeat; i > 0; --i)
            *dst++ += (b - (*src1++ + *src2++)) >> e;
      }
      else { // general case
        if (synthesis)
          for (ui32 i = repeat; i > 0; --i)
            *dst++ -= (b + a * (*src1++ + *src2++)) >> e;
        else
          for (ui32 i = repeat; i > 0; --i)
            *dst++ += (b + a * (*src1++ + *src2++)) >> e;
      }
    }

    static
    void gen_rev_vert_step64(const lifting_step* s, const line_buf* sig, 
                             const line_buf* other, const line_buf* aug, 
                             ui32 repeat, bool synthesis)
    {
      const si64 a = s->rev.Aatk;
      const si64 b = s->rev.Batk;
      const ui8 e = s->rev.Eatk;
 
      si64* dst = aug->i64;
      const si64* src1 = sig->i64, * src2 = other->i64;
      // The general definition of the wavelet in Part 2 is slightly 
      // different to part 2, although they are mathematically equivalent
      // here, we identify the simpler form from Part 1 and employ them
      if (a == 1)
      { // 5/3 update and any case with a == 1
        if (synthesis)
          for (ui32 i = repeat; i > 0; --i)
            *dst++ -= (b + *src1++ + *src2++) >> e;
        else
          for (ui32 i = repeat; i > 0; --i)
            *dst++ += (b + *src1++ + *src2++) >> e;
      }
      else if (a == -1 && b == 1 && e == 1)
      { // 5/3 predict
        if (synthesis)
          for (ui32 i = repeat; i > 0; --i)
            *dst++ += (*src1++ + *src2++) >> e;
        else
          for (ui32 i = repeat; i > 0; --i)
            *dst++ -= (*src1++ + *src2++) >> e;
      }
      else if (a == -1)
      { // any case with a == -1, which is not 5/3 predict
        if (synthesis)
          for (ui32 i = repeat; i > 0; --i)
            *dst++ -= (b - (*src1++ + *src2++)) >> e;
        else
          for (ui32 i = repeat; i > 0; --i)
            *dst++ += (b - (*src1++ + *src2++)) >> e;
      }
      else { // general case
        if (synthesis)
          for (ui32 i = repeat; i > 0; --i)
            *dst++ -= (b + a * (*src1++ + *src2++)) >> e;
        else
          for (ui32 i = repeat; i > 0; --i)
            *dst++ += (b + a * (*src1++ + *src2++)) >> e;
      }
    }

    void gen_rev_vert_step(const lifting_step* s, const line_buf* sig, 
                           const line_buf* other, const line_buf* aug, 
                           ui32 repeat, bool synthesis)
    {
      if (((sig != NULL) && (sig->flags & line_buf::LFT_32BIT)) || 
          ((aug != NULL) && (aug->flags & line_buf::LFT_32BIT)) ||
          ((other != NULL) && (other->flags & line_buf::LFT_32BIT))) 
      {
        assert((sig == NULL || sig->flags & line_buf::LFT_32BIT) &&
               (other == NULL || other->flags & line_buf::LFT_32BIT) && 
               (aug == NULL || aug->flags & line_buf::LFT_32BIT));
        gen_rev_vert_step32(s, sig, other, aug, repeat, synthesis);
      }
      else 
      {
        assert((sig == NULL || sig->flags & line_buf::LFT_64BIT) &&
               (other == NULL || other->flags & line_buf::LFT_64BIT) && 
               (aug == NULL || aug->flags & line_buf::LFT_64BIT));
        gen_rev_vert_step64(s, sig, other, aug, repeat, synthesis);
      }
    }

    static
    void gen_rev_horz_ana32(const param_atk* atk, const line_buf* ldst, 
                            const line_buf* hdst, const line_buf* src, 
                            ui32 width, bool even)
    {
      if (width > 1)
      {
        // combine both lsrc and hsrc into dst
        si32* dph = hdst->i32;
        si32* dpl = ldst->i32;
        si32* sp = src->i32;
        ui32 w = width;
        if (!even)
        {
          *dph++ = *sp++; --w;
        }
        for (; w > 1; w -= 2)
        {
          *dpl++ = *sp++; *dph++ = *sp++;
        }
        if (w)
        {
          *dpl++ = *sp++; --w;
        }
 
        si32* hp = hdst->i32, * lp = ldst->i32;
        ui32 l_width = (width + (even ? 1 : 0)) >> 1;  // low pass
        ui32 h_width = (width + (even ? 0 : 1)) >> 1;  // high pass
        ui32 num_steps = atk->get_num_steps();
        for (ui32 j = num_steps; j > 0; --j)
        {
          // first lifting step
          const lifting_step* s = atk->get_step(j - 1);
          const si32 a = s->rev.Aatk;
          const si32 b = s->rev.Batk;
          const ui8 e = s->rev.