ojph_params.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_params.cpp
// Author: Aous Naman
// Date: 28 August 2019
//***************************************************************************/
 
#define _USE_MATH_DEFINES
#include <cmath>
#include <climits>
#include <new>
 
#include "ojph_arch.h"
#include "ojph_base.h"
#include "ojph_file.h"
#include "ojph_params.h"
 
#include "ojph_params_local.h"
#include "ojph_message.h"
 
namespace ojph {

  //
  //
  //
  //
  //

  void param_siz::set_image_extent(point dims)
  {
    state->set_image_extent(dims);
  }

  void param_siz::set_tile_size(size s)
  {
    state->set_tile_size(s);
  }

  void param_siz::set_image_offset(point offset)
  {
    state->set_image_offset(offset);
  }

  void param_siz::set_tile_offset(point offset)
  {
    state->set_tile_offset(offset);
  }

  void param_siz::set_num_components(ui32 num_comps)
  {
    state->set_num_components(num_comps);
  }

  void param_siz::set_component(ui32 comp_num, const point& downsampling,
                                ui32 bit_depth, bool is_signed)
  {
    state->set_comp_info(comp_num, downsampling, bit_depth, is_signed);
  }

  point param_siz::get_image_extent() const
  {
    return point(state->Xsiz, state->Ysiz);
  }

  point param_siz::get_image_offset() const
  {
    return point(state->XOsiz, state->YOsiz);
  }

  size param_siz::get_tile_size() const
  {
    return size(state->XTsiz, state->YTsiz);
  }

  point param_siz::get_tile_offset() const
  {
    return point(state->XTOsiz, state->YTOsiz);
  }

  ui32 param_siz::get_num_components() const
  {
    return state->Csiz;
  }

  ui32 param_siz::get_bit_depth(ui32 comp_num) const
  {
    return state->get_bit_depth(comp_num);
  }

  bool param_siz::is_signed(ui32 comp_num) const
  {
    return state->is_signed(comp_num);
  }

  point param_siz::get_downsampling(ui32 comp_num) const
  {
    return state->get_downsampling(comp_num);
  }

  ui32 param_siz::get_recon_width(ui32 comp_num) const
  {
    return state->get_recon_width(comp_num);
  }

  ui32 param_siz::get_recon_height(ui32 comp_num) const
  {
    return state->get_recon_height(comp_num);
  }

  //
  //
  //
  //
  //

  void param_cod::set_num_decomposition(ui32 num_decompositions)
  {
    if (num_decompositions > 32)
      OJPH_ERROR(0x00050001,
        "maximum number of decompositions cannot exceed 32");
    state->SPcod.num_decomp = (ui8)num_decompositions;
  }

  void param_cod::set_block_dims(ui32 width, ui32 height)
  {
    ui32 log_width = 31 - count_leading_zeros(width);
    ui32 log_height = 31 - count_leading_zeros(height);
    if (width == 0 || width != (1u << log_width)
      || height == 0 || height != (1u << log_height)
      || log_width < 2 || log_height < 2
      || log_width + log_height > 12)
      OJPH_ERROR(0x00050011, "incorrect code block dimensions");
    state->SPcod.block_width = (ui8)(log_width - 2);
    state->SPcod.block_height = (ui8)(log_height - 2);
  }

  void param_cod::set_precinct_size(int num_levels, size* precinct_size)
  {
    if (num_levels == 0 || precinct_size == NULL)
      state->Scod &= 0xFE;
    else
    {
      state->Scod |= 1;
      for (int i = 0; i <= state->SPcod.num_decomp; ++i)
      {
        size t = precinct_size[i < num_levels ? i : num_levels - 1];
 
        ui32 PPx = 31 - count_leading_zeros(t.w);
        ui32 PPy = 31 - count_leading_zeros(t.h);
        if (t.w == 0 || t.h == 0)
          OJPH_ERROR(0x00050021, "precinct width or height cannot be 0");
        if (t.w != (1u<<PPx) || t.h != (1u<<PPy))
          OJPH_ERROR(0x00050022,
            "precinct width and height should be a power of 2");
        if (PPx > 15 || PPy > 15)
          OJPH_ERROR(0x00050023, "precinct size is too large");
        if (i > 0 && (PPx == 0 || PPy == 0))
          OJPH_ERROR(0x00050024, "precinct size is too small");
        state->SPcod.precinct_size[i] = (ui8)(PPx | (PPy << 4));
      }
    }
  }

  void param_cod::set_progression_order(const char *name)
  {
    int prog_order = 0;
    size_t len = strlen(name);
    if (len == 4)
    {
      if (strncmp(name, OJPH_PO_STRING_LRCP, 4) == 0)
        prog_order = OJPH_PO_LRCP;
      else if (strncmp(name, OJPH_PO_STRING_RLCP, 4) == 0)
        prog_order = OJPH_PO_RLCP;
      else if (strncmp(name, OJPH_PO_STRING_RPCL, 4) == 0)
        prog_order = OJPH_PO_RPCL;
      else if (strncmp(name, OJPH_PO_STRING_PCRL, 4) == 0)
        prog_order = OJPH_PO_PCRL;
      else if (strncmp(name, OJPH_PO_STRING_CPRL, 4) == 0)
        prog_order = OJPH_PO_CPRL;
      else
        OJPH_ERROR(0x00050031, "unknown progression order");
    }
    else
      OJPH_ERROR(0x00050032, "improper progression order");
 
 
    state->SGCod.prog_order = (ui8)prog_order;
  }

  void param_cod::set_color_transform(bool color_transform)
  {
    state->employ_color_transform(color_transform ? 1 : 0);
  }

  void param_cod::set_reversible(bool reversible)
  {
    state->set_reversible(reversible);
  }

  void param_cod::set_num_decomposition(ui32 comp_idx, ui32 num_decompositions)
  {
    local::param_cod* cdp = state->get_or_add_coc(comp_idx);
    ojph::param_cod(cdp).set_num_decomposition(num_decompositions);
  }

  void param_cod::set_block_dims(ui32 comp_idx, ui32 width, ui32 height)
  {
    local::param_cod* cdp = state->get_or_add_coc(comp_idx);
    ojph::param_cod(cdp).set_block_dims(width, height);
  }

  void param_cod::set_precinct_size(ui32 comp_idx, int num_levels,
                                    size* precinct_size)
  {
    local::param_cod* cdp = state->get_or_add_coc(comp_idx);
    ojph::param_cod(cdp).set_precinct_size(num_levels, precinct_size);
  }

  void param_cod::set_reversible(ui32 comp_idx, bool reversible)
  {
    local::param_cod* cdp = state->get_or_add_coc(comp_idx);
    ojph::param_cod(cdp).set_reversible(reversible);
  }

  ui32 param_cod::get_num_decompositions() const
  {
    return state->get_num_decompositions();
  }

  size param_cod::get_block_dims() const
  {
    return state->get_block_dims();
  }

  size param_cod::get_log_block_dims() const
  {
    return state->get_log_block_dims();
  }

  bool param_cod::is_reversible() const
  {
    return state->is_reversible();
  }

  size param_cod::get_precinct_size(ui32 level_num) const
  {
    return state->get_precinct_size(level_num);
  }

  size param_cod::get_log_precinct_size(ui32 level_num) const
  {
    return state->get_log_precinct_size(level_num);
  }

  int param_cod::get_progression_order() const
  {
    return state->SGCod.prog_order;
  }

  const char* param_cod::get_progression_order_as_string() const
  {
    if (state->SGCod.prog_order == OJPH_PO_LRCP)
      return OJPH_PO_STRING_LRCP;
    else if (state->SGCod.prog_order == OJPH_PO_RLCP)
      return OJPH_PO_STRING_RLCP;
    else if (state->SGCod.prog_order == OJPH_PO_RPCL)
      return OJPH_PO_STRING_RPCL;
    else if (state->SGCod.prog_order == OJPH_PO_PCRL)
      return OJPH_PO_STRING_PCRL;
    else if (state->SGCod.prog_order == OJPH_PO_CPRL)
      return OJPH_PO_STRING_CPRL;
    else
      assert(0);
    return "";
  }

  int param_cod::get_num_layers() const
  {
    return state->SGCod.num_layers;
  }

  bool param_cod::is_using_color_transform() const
  {
    return state->is_employing_color_transform();
  }

  bool param_cod::packets_may_use_sop() const
  {
    return state->packets_may_use_sop();
  }

  bool param_cod::packets_use_eph() const
  {
    return state->packets_use_eph();
  }

  bool param_cod::get_block_vertical_causality() const
  {
    return state->get_block_vertical_causality();
  }

  ui32 param_cod::get_num_decompositions(ui32 comp_idx) const
  { return state->get_coc(comp_idx)->get_num_decompositions(); }

  size param_cod::get_block_dims(ui32 comp_idx) const
  { return state->get_coc(comp_idx)->get_block_dims(); }

  size param_cod::get_log_block_dims(ui32 comp_idx) const
  { return state->get_coc(comp_idx)->get_log_block_dims(); }

  bool param_cod::is_reversible(ui32 comp_idx) const
  { return state->get_coc(comp_idx)->is_reversible(); }

  size param_cod::get_precinct_size(ui32 comp_idx, ui32 level_num) const
  { return state->get_coc(comp_idx)->get_precinct_size(level_num); }

  size param_cod::get_log_precinct_size(ui32 comp_idx, ui32 level_num) const
  { return state->get_coc(comp_idx)->get_log_precinct_size(level_num); }

  bool param_cod::get_block_vertical_causality(ui32 comp_idx) const
  { return state->get_coc(comp_idx)->get_block_vertical_causality(); }
 

  //
  //
  //
  //
  //

  void param_qcd::set_irrev_quant(float delta)
  {
    state->set_delta(delta);
  }

  void param_qcd::set_irrev_quant(float delta, size_t num_weights,
                                  const float weights[])
  {
    state->set_delta(delta, num_weights, weights);
  }

  void param_qcd::set_qfactor(float qfactor) {
    state->set_qfactor(qfactor);
  }

  void param_qcd::set_irrev_quant(ui32 comp_idx, float delta)
  {
    state->set_delta(comp_idx, delta);
  }

  void param_qcd::set_irrev_quant(ui32 comp_idx, float delta,
                                  size_t num_weights, const float weights[])
  {
    state->set_delta(comp_idx, delta, num_weights, weights);
  }

  void param_qcd::set_qfactor(ui32 comp_idx, comp_type ctype, float qfactor)
  {
    state->set_qfactor(comp_idx, ctype, qfactor);
  }

  //
  //
  //
  //
  //

  void param_nlt::set_nonlinear_transform(ui32 comp_num, ui8 nl_type)
  {
    state->set_nonlinear_transform(comp_num, nl_type);
  }

  void param_nlt::set_nonlinear_transform(ui32 comp_num,
                                          ui8 decoded_bit_depth,
                                          bool decoded_signedness,
                                          ui32 d_min, ui32 d_max, ui8 pt_val,
                                          ui16 num_points, void* points,
                                          ui8 nl_type,
                                          bool approx_encoding_nlt,
                                          ui32 num_enc_points)
  {
    state->set_nonlinear_transform(comp_num, decoded_bit_depth,
      decoded_signedness, d_min, d_max, pt_val, num_points, points, nl_type,
      approx_encoding_nlt, num_enc_points);
  }

  bool param_nlt::get_nonlinear_transform(ui32 comp_num,
                                          ui8& decoded_bit_depth,
                                          bool& decoded_signedness,
                                          ui8& nl_type) const
  {
    return state->get_nonlinear_transform(comp_num, decoded_bit_depth,
      decoded_signedness, nl_type);
  }

  //
  //
  //
  //
  //

  void comment_exchange::set_string(const char* str)
  {
    size_t t = strlen(str);
    if (len > 65531)
      OJPH_ERROR(0x000500C1,
        "COM marker string length cannot be larger than 65531");
    this->data = str;
    this->len = (ui16)t;
    this->Rcom = 1;
  }

  void comment_exchange::set_data(const char* data, ui16 len)
  {
    if (len > 65531)
      OJPH_ERROR(0x000500C2,
        "COM marker string length cannot be larger than 65531");
    this->data = data;
    this->len = len;
    this->Rcom = 0;
  }

  //
  //
  //                                LOCAL
  //
  //
 
  namespace local {

    //static
    class sqrt_energy_gains
    {
    public:
      static float get_gain_l(ui32 num_decomp, bool reversible)
      { return reversible ? gain_5x3_l[num_decomp] : gain_9x7_l[num_decomp]; }
      static float get_gain_h(ui32 num_decomp, bool reversible)
      { return reversible ? gain_5x3_h[num_decomp] : gain_9x7_h[num_decomp]; }
 
    private:
      static const float gain_9x7_l[34];
      static const float gain_9x7_h[34];
      static const float gain_5x3_l[34];
      static const float gain_5x3_h[34];
    };

    const float sqrt_energy_gains::gain_9x7_l[34] = { 1.0000e+00f,
      1.4021e+00f, 2.0304e+00f, 2.9012e+00f, 4.1153e+00f, 5.8245e+00f,
      8.2388e+00f, 1.1652e+01f, 1.6479e+01f, 2.3304e+01f, 3.2957e+01f,
      4.6609e+01f, 6.5915e+01f, 9.3217e+01f, 1.3183e+02f, 1.8643e+02f,
      2.6366e+02f, 3.7287e+02f, 5.2732e+02f, 7.4574e+02f, 1.0546e+03f,
      1.4915e+03f, 2.1093e+03f, 2.9830e+03f, 4.2185e+03f, 5.9659e+03f,
      8.4371e+03f, 1.1932e+04f, 1.6874e+04f, 2.3864e+04f, 3.3748e+04f,
      4.7727e+04f, 6.7496e+04f, 9.5454e+04f };
    const float sqrt_energy_gains::gain_9x7_h[34] = { 1.4425e+00f,
      1.9669e+00f, 2.8839e+00f, 4.1475e+00f, 5.8946e+00f, 8.3472e+00f,
      1.1809e+01f, 1.6701e+01f, 2.3620e+01f, 3.3403e+01f, 4.7240e+01f,
      6.6807e+01f, 9.4479e+01f, 1.3361e+02f, 1.8896e+02f, 2.6723e+02f,
      3.7792e+02f, 5.3446e+02f, 7.5583e+02f, 1.0689e+03f, 1.5117e+03f,
      2.1378e+03f, 3.0233e+03f, 4.2756e+03f, 6.0467e+03f, 8.5513e+03f,
      1.2093e+04f, 1.7103e+04f, 2.4187e+04f, 3.4205e+04f, 4.8373e+04f,
      6.8410e+04f, 9.6747e+04f, 1.3682e+05f };
    const float sqrt_energy_gains::gain_5x3_l[34] = { 1.0000e+00f,
      1.2247e+00f, 1.3229e+00f, 1.5411e+00f, 1.7139e+00f, 1.9605e+00f,
      2.2044e+00f, 2.5047e+00f, 2.8277e+00f, 3.2049e+00f, 3.6238e+00f,
      4.1033e+00f, 4.6423e+00f, 5.2548e+00f, 5.9462e+00f, 6.7299e+00f,
      7.6159e+00f, 8.6193e+00f, 9.7544e+00f, 1.1039e+01f, 1.2493e+01f,
      1.4139e+01f, 1.6001e+01f, 1.8108e+01f, 2.0493e+01f, 2.3192e+01f,
      2.6246e+01f, 2.9702e+01f, 3.3614e+01f, 3.8041e+01f, 4.3051e+01f,
      4.8721e+01f, 5.5138e+01f, 6.2399e+01f };
    const float sqrt_energy_gains::gain_5x3_h[34] = { 1.0458e+00f,
      1.3975e+00f, 1.4389e+00f, 1.7287e+00f, 1.8880e+00f, 2.1841e+00f,
      2.4392e+00f, 2.7830e+00f, 3.1341e+00f, 3.5576e+00f, 4.0188e+00f,
      4.5532e+00f, 5.1494e+00f, 5.8301e+00f, 6.5963e+00f, 7.4663e+00f,
      8.4489e+00f, 9.5623e+00f, 1.0821e+01f, 1.2247e+01f, 1.3860e+01f,
      1.5685e+01f, 1.7751e+01f, 2.0089e+01f, 2.2735e+01f, 2.5729e+01f,
      2.9117e+01f, 3.2952e+01f, 3.7292e+01f, 4.2203e+01f, 4.7761e+01f,
      5.4051e+01f, 6.1170e+01f, 6.9226e+01f };

    //static
    class bibo_gains
    {
    public:
      static float get_bibo_gain_l(ui32 num_decomp, bool reversible)
      { return reversible ? gain_5x3_l[num_decomp] : gain_9x7_l[num_decomp]; }
      static float get_bibo_gain_h(ui32 num_decomp, bool reversible)
      { return reversible ? gain_5x3_h[num_decomp] : gain_9x7_h[num_decomp]; }
 
    private:
      static const float gain_9x7_l[34];
      static const float gain_9x7_h[34];
      static const float gain_5x3_l[34];
      static const float gain_5x3_h[34];
    };

    const float bibo_gains::gain_9x7_l[34] = { 1.0000e+00f, 1.3803e+00f,
      1.3328e+00f, 1.3067e+00f, 1.3028e+00f, 1.3001e+00f, 1.2993e+00f,
      1.2992e+00f, 1.2992e+00f, 1.2992e+00f, 1.2992e+00f, 1.2992e+00f,
      1.2992e+00f, 1.2992e+00f, 1.2992e+00f, 1.2992e+00f, 1.2992e+00f,
      1.2992e+00f, 1.2992e+00f, 1.2992e+00f, 1.2992e+00f, 1.2992e+00f,
      1.2992e+00f, 1.2992e+00f, 1.2992e+00f, 1.2992e+00f, 1.2992e+00f,
      1.2992e+00f, 1.2992e+00f, 1.2992e+00f, 1.2992e+00f, 1.2992e+00f,
      1.2992e+00f, 1.2992e+00f };
    const float bibo_gains::gain_9x7_h[34] = { 1.2976e+00f, 1.3126e+00f,
      1.2757e+00f, 1.2352e+00f, 1.2312e+00f, 1.2285e+00f, 1.2280e+00f,
      1.2278e+00f, 1.2278e+00f, 1.2278e+00f, 1.2278e+00f, 1.2278e+00f,
      1.2278e+00f, 1.2278e+00f, 1.2278e+00f, 1.2278e+00f, 1.2278e+00f,
      1.2278e+00f, 1.2278e+00f, 1.2278e+00f, 1.2278e+00f, 1.2278e+00f,
      1.2278e+00f, 1.2278e+00f, 1.2278e+00f, 1.2278e+00f, 1.2278e+00f,
      1.2278e+00f, 1.2278e+00f, 1.2278e+00f, 1.2278e+00f, 1.2278e+00f,
      1.2278e+00f, 1.2278e+00f };
    const float bibo_gains::gain_5x3_l[34] = { 1.0000e+00f, 1.5000e+00f,
      1.6250e+00f, 1.6875e+00f, 1.6963e+00f, 1.7067e+00f, 1.7116e+00f,
      1.7129e+00f, 1.7141e+00f, 1.7145e+00f, 1.7151e+00f, 1.7152e+00f,
      1.7155e+00f, 1.7155e+00f, 1.7156e+00f, 1.7156e+00f, 1.7156e+00f,
      1.7156e+00f, 1.7156e+00f, 1.7156e+00f, 1.7156e+00f, 1.7156e+00f,
      1.7156e+00f, 1.7156e+00f, 1.7156e+00f, 1.7156e+00f, 1.7156e+00f,
      1.7156e+00f, 1.7156e+00f, 1.7156e+00f, 1.7156e+00f, 1.7156e+00f,
      1.7156e+00f, 1.7156e+00f };
    const float bibo_gains::gain_5x3_h[34] = { 2.0000e+00f, 2.5000e+00f,
      2.7500e+00f, 2.8047e+00f, 2.8198e+00f, 2.8410e+00f, 2.8558e+00f,
      2.8601e+00f, 2.8628e+00f, 2.8656e+00f, 2.8662e+00f, 2.8667e+00f,
      2.8669e+00f, 2.8670e+00f, 2.8671e+00f, 2.8671e+00f, 2.8671e+00f,
      2.8671e+00f, 2.8671e+00f, 2.8671e+00f, 2.8671e+00f, 2.8671e+00f,
      2.8671e+00f, 2.8671e+00f, 2.8671e+00f, 2.8671e+00f, 2.8671e+00f,
      2.8671e+00f, 2.8671e+00f, 2.8671e+00f, 2.8671e+00f, 2.8671e+00f,
      2.8671e+00f, 2.8671e+00f };

    //static
    class visual_weights
    {
      public:
        enum colour_format : ui32
        {
          VW_COLOUR_FORMAT_400 = 0,
          VW_COLOUR_FORMAT_420 = 1,
          VW_COLOUR_FORMAT_422 = 2,
          VW_COLOUR_FORMAT_444 = 3,
          VW_COLOUR_FORMAT_ERROR = 4
        };
        using comp_type = param_qcd::comp_type;
 
    public:
      static colour_format get_format(const point& component_subsampling)
      {
        if (component_subsampling.x == 2 && component_subsampling.y == 2)
          return VW_COLOUR_FORMAT_420;
        else if (component_subsampling.x == 2 && component_subsampling.y == 1)
          return VW_COLOUR_FORMAT_422;
        else if (component_subsampling.x == 1 && component_subsampling.y == 1)
          return VW_COLOUR_FORMAT_444;
        else
          return VW_COLOUR_FORMAT_ERROR;
      }
 
    public:
      static const float* get_weights(ui32 format, ui32 comp_type)
      {
        if (comp_type == comp_type::OJPH_COMP_Y)
          return y;
        else if (comp_type == comp_type::OJPH_COMP_CB)
        {
          if (format == VW_COLOUR_FORMAT_420)        return cb420;
          else if (format == VW_COLOUR_FORMAT_422)   return cb422;
          else if (format == VW_COLOUR_FORMAT_444)   return cb444;
          else {
            assert(0);
            return y;
          }
        }
        else if (comp_type == comp_type::OJPH_COMP_CR)
        {
          if (format == VW_COLOUR_FORMAT_420)        return cr420;
          else if (format == VW_COLOUR_FORMAT_422)   return cr422;
          else if (format == VW_COLOUR_FORMAT_444)   return cr444;
          else {
            assert(0);
            return y;
          }
        }
        else {
          assert(0);
          return y;
        }
      }
 
      static const float* get_no_weights()
      { return no_weights; }
 
      static float get_weight(const float *v, ui32 decomposition_level,
                              ui32 subband_idx)
 
      {
        if (subband_idx == 0)
          return v[18];
        else {
          assert(subband_idx >= 1 && subband_idx <= 3);
          assert(decomposition_level > 0);
          decomposition_level = ojph_min(decomposition_level, 6);
          ui32 index = (decomposition_level - 1) * 3 + (3 - subband_idx);
          return v[index];
        }
      }

      static float get_gain(ui32 comp_type)
      {
        if (comp_type == comp_type::OJPH_COMP_Y)
          return 1.0f;
        else if (comp_type == comp_type::OJPH_COMP_CB)
          return 1.8051f / 1.7321f;
        else if (comp_type == comp_type::OJPH_COMP_CR)
          return 1.5734f / 1.7321f;
        else {
          assert(0);
          return 0.0f;
        }
      }

      static float get_delta_ref(float qfactor, ui32 bit_depth,
                                 float& power)
      {
        // returns delta_ref & power to be used with visual weights
        constexpr float t0       = 65.0f;
        constexpr float t1       = 97.0f;
        constexpr float alpha_t0 = 0.04f, alpha_t1 = 0.10f;
        constexpr float m_t0     = 2.0f * (1.0f - t0 / 100.0f);
        constexpr float m_t1     = 2.0f * (1.0f - t1 / 100.0f);
 
        float m_q;
        if (qfactor < 50.0f)
          m_q = 50.0f / qfactor;
        else
          m_q = 2.0f * (1.0f - qfactor / 100.0f);
 
        float alpha_q;
        if (qfactor <= t0)
        {
          power = 1.0f;
          alpha_q = alpha_t0;
        }
        else if (qfactor < t1)
        {
          power = std::log(m_q) - std::log(m_t1);
          power /= std::log(m_t0) - std::log(m_t1);
          alpha_q = alpha_t1 * std::pow(alpha_t0 / alpha_t1, power);
        }
        else
        {
          power = 0.0f;
          alpha_q = alpha_t1;
        }
        const float eps = std::sqrt(0.5f) * std::ldexp(1.0f, -(int)bit_depth);
        return alpha_q * m_q + eps;
      }
 
    private:
      static const float cb420[19];
      static const float cr420[19];
      static const float cb422[19];
      static const float cr422[19];
      static const float cb444[19];
      static const float cr444[19];
      static const float y[19];
      static const float no_weights[19];
    };

    const float visual_weights::cb420[19] = {
      0.2724f, 0.5128f, 0.5128f,              // level 1
      0.6692f, 0.9382f, 0.9382f,              // level 2
      1.0888f, 1.3046f, 1.3046f,              // level 3
      1.4156f, 1.5594f, 1.5594f,              // level 4
      2.0f,    2.0f,    2.0f,                 // level 5
      2.0f,    2.0f,    2.0f,    2.0f};       // level 6 + LL
    const float visual_weights::cr420[19] = {
      0.5196f, 0.8260f, 0.8260f,              // level 1
      1.0080f, 1.2928f, 1.2928f,              // level 2
      1.4440f, 1.6508f, 1.6508f,              // level 3
      1.7538f, 1.8848f, 1.8848f,              // level 4
      2.0f,    2.0f,    2.0f,                 // level 5
      2.0f,    2.0f,    2.0f,    2.0f};       // level 6 + LL
    const float visual_weights::cb422[19] = {
      0.1220f, 0.1220f, 0.3626f,              // level 1
      0.3626f, 0.3626f, 0.6634f,              // level 2
      0.6634f, 0.6634f, 0.9225f,              // level 3
      0.9225f, 0.9225f, 1.1027f,              // level 4
      1.1027f, 1.1027f, 1.4142f,              // level 5
      1.4142f, 1.4142f, 1.4142f, 1.4142f};    // level 6 + LL
    const float visual_weights::cr422[19] = {
      0.2595f, 0.2595f, 0.5841f,              // level 1
      0.5841f, 0.5841f, 0.9141f,              // level 2
      0.9141f, 0.9141f, 1.1673f,              // level 3
      1.1673f, 1.1673f, 1.3328f,              // level 4
      1.3328f, 1.3328f, 1.4142f,              // level 5
      1.4142f, 1.4142f, 1.4142f, 1.4142f};    // level 6 + LL
    const float visual_weights::cb444[19] = {
      0.0263f, 0.0863f, 0.0863f,              // level 1
      0.1362f, 0.2564f, 0.2564f,              // level 2
      0.3346f, 0.4691f, 0.4691f,              // level 3
      0.5444f, 0.6523f, 0.6523f,              // level 4
      0.7078f, 0.7797f, 0.7797f,              // level 5
      1.0f,    1.0f,    1.0f,    1.0f};       // level 6 + LL
    const float visual_weights::cr444[19] = {
      0.0773f, 0.1835f, 0.1835f,              // level 1
      0.2598f, 0.4130f, 0.4130f,              // level 2
      0.5040f, 0.6464f, 0.6464f,              // level 3
      0.7220f, 0.8254f, 0.8254f,              // level 4
      0.8769f, 0.9424f, 0.9424f,              // level 5
      1.0f,    1.0f,    1.0f,    1.0f};       // level 6 + LL
    const float visual_weights::y[19]     = {
      0.0901f, 0.2758f, 0.2758f,              // level 1
      0.7018f, 0.8378f, 0.8378f,              // level 2
      1.0f,    1.0f,    1.0f,                 // level 3
      1.0f,    1.0f,    1.0f,                 // level 4
      1.0f,    1.0f,    1.0f,                 // level 5
      1.0f,    1.0f,    1.0f,    1.0f};       // level 6 + LL
    const float visual_weights::no_weights[19] = {
      1.0f,    1.0f,    1.0f,                 // level 1
      1.0f,    1.0f,    1.0f,                 // level 2
      1.0f,    1.0f,    1.0f,                 // level 3
      1.0f,    1.0f,    1.0f,                 // level 4
      1.0f,    1.0f,    1.0f,                 // level 5
      1.0f,    1.0f,    1.0f,    1.0f};       // level 6 + LL
 

    //
    //
    //
    //
    //

    bool param_siz::write(outfile_base *file)
    {
      //marker size excluding header
      Lsiz = (ui16)(38 + 3 * Csiz);
 
      ui8  buf1;
      ui16 buf2;
      ui32 buf4;
      bool result = true;
 
      buf2 = JP2K_MARKER::SIZ;
      buf2 = swap_bytes_if_le(buf2);
      result &= file->write(&buf2, sizeof(ui16)) == sizeof(ui16);
      buf2 = swap_bytes_if_le(Lsiz);
      result &= file->write(&buf2, sizeof(ui16)) == sizeof(ui16);
      buf2 = swap_bytes_if_le(Rsiz);
      result &= file->write(&buf2, sizeof(ui16)) == sizeof(ui16);
      buf4 = swap_bytes_if_le(Xsiz);
      result &= file->write(&buf4, sizeof(ui32)) == sizeof(ui32);
      buf4 = swap_bytes_if_le(Ysiz);
      result &= file->write(&buf4, sizeof(ui32)) == sizeof(ui32);
      buf4 = swap_bytes_if_le(XOsiz);
      result &= file->write(&buf4, sizeof(ui32)) == sizeof(ui32);
      buf4 = swap_bytes_if_le(YOsiz);
      result &= file->write(&buf4, sizeof(ui32)) == sizeof(ui32);
      buf4 = swap_bytes_if_le(XTsiz);
      result &= file->write(&buf4, sizeof(ui32)) == sizeof(ui32);
      buf4 = swap_bytes_if_le(YTsiz);
      result &= file->write(&buf4, sizeof(ui32)) == sizeof(ui32);
      buf4 = swap_bytes_if_le(XTOsiz);
      result &= file->write(&buf4, sizeof(ui32)) == sizeof(ui32);
      buf4 = swap_bytes_if_le(YTOsiz);
      result &= file->write(&buf4, sizeof(ui32)) == sizeof(ui32);
      buf2 = swap_bytes_if_le(Csiz);
      result &= file->write(&buf2, sizeof(ui16)) == sizeof(ui16);
      for (int c = 0; c < Csiz; ++c)
      {
        buf1 = cptr[c].SSiz;
        result &= file->write(&buf1, sizeof(ui8)) == sizeof(ui8);
        buf1 = cptr[c].XRsiz;
        result &= file->write(&buf1, sizeof(ui8)) == sizeof(ui8);
        buf1 = cptr[c].YRsiz;
        result &= file->write(&buf1, sizeof(ui8)) == sizeof(ui8);
      }
 
      return result;
    }

    void param_siz::read(infile_base *file)
    {
      if (file->read(&Lsiz, 2) != 2)
        OJPH_ERROR(0x00050041, "error reading SIZ marker");
      Lsiz = swap_bytes_if_le(Lsiz);
      int num_comps = (Lsiz - 38) / 3;
      if (Lsiz != 38 + 3 * num_comps)
        OJPH_ERROR(0x00050042, "error in SIZ marker length");
      if (file->read(&Rsiz, 2) != 2)
        OJPH_ERROR(0x00050043, "error reading SIZ marker");
      Rsiz = swap_bytes_if_le(Rsiz);
      if ((Rsiz & 0x4000) == 0)
        OJPH_ERROR(0x00050044,
          "Rsiz bit 14 is not set (this is not a JPH file)");
      if ((Rsiz & 0x8000) != 0 && (Rsiz & 0xD5F) != 0)
        OJPH_WARN(0x00050001, "Rsiz in SIZ has unimplemented fields");
      if (file->read(&Xsiz, 4) != 4)
        OJPH_ERROR(0x00050045, "error reading SIZ marker");
      Xsiz = swap_bytes_if_le(Xsiz);
      if (file->read(&Ysiz, 4) != 4)
        OJPH_ERROR(0x00050046, "error reading SIZ marker");
      Ysiz = swap_bytes_if_le(Ysiz);
      ui32 t_XOsiz, t_YOsiz;
      if (file->read(&t_XOsiz, 4) != 4)
        OJPH_ERROR(0x00050047, "error reading SIZ marker");
      if (file->read(&t_YOsiz, 4) != 4)
        OJPH_ERROR(0x00050048, "error reading SIZ marker");
      set_image_offset(point(
        swap_bytes_if_le(t_XOsiz),
        swap_bytes_if_le(t_YOsiz)));
      ui32 t_XTsiz, t_YTsiz;
      if (file->read(&t_XTsiz, 4) != 4)
        OJPH_ERROR(0x00050049, "error reading SIZ marker");
      if (file->read(&t_YTsiz, 4) != 4)
        OJPH_ERROR(0x0005004A, "error reading SIZ marker");
      set_tile_size(size(
        swap_bytes_if_le(t_XTsiz),
        swap_bytes_if_le(t_YTsiz)));
      ui32 t_XTOsiz, t_YTOsiz;
      if (file->read(&t_XTOsiz, 4) != 4)
        OJPH_ERROR(0x0005004B, "error reading SIZ marker");
      if (file->read(&t_YTOsiz, 4) != 4)
        OJPH_ERROR(0x0005004C, "error reading SIZ marker");
      set_tile_offset(point(
        swap_bytes_if_le(t_XTOsiz),
        swap_bytes_if_le(t_YTOsiz)));
      if (file->read(&Csiz, 2) != 2)
        OJPH_ERROR(0x0005004D, "error reading SIZ marker");
      Csiz = swap_bytes_if_le(Csiz);
      if (Csiz != num_comps)
        OJPH_ERROR(0x0005004E, "Csiz does not match the SIZ marker size");
      if (Csiz == 0)
        OJPH_ERROR(0x0005004F, "Wrong Csiz value of 0 in SIZ marker segment");
      set_num_components(Csiz);
      for (int c = 0; c < Csiz; ++c)
      {
        if (file->read(&cptr[c].SSiz, 1) != 1)
          OJPH_ERROR(0x00050051, "error reading SIZ marker");
        if (file->read(&cptr[c].XRsiz, 1) != 1)
          OJPH_ERROR(0x00050052, "error reading SIZ marker");
        if (file->read(&cptr[c].YRsiz, 1) != 1)
          OJPH_ERROR(0x00050053, "error reading SIZ marker");
        if ((cptr[c].SSiz & 0x7F) > 37)
          OJPH_ERROR(0x00050054, "Wrong SIZ-SSiz value of %d", cptr[c].SSiz);
        if (cptr[c].XRsiz == 0)
          OJPH_ERROR(0x00050055, "Wrong SIZ-XRsiz value of %d", cptr[c].XRsiz);
        if (cptr[c].YRsiz == 0)
          OJPH_ERROR(0x00050056, "Wrong SIZ-YRsiz value of %d", cptr[c].YRsiz);
      }
 
      ws_kern_support_needed = (Rsiz & 0x20) != 0;
      dfs_support_needed = (Rsiz & 0x80) != 0;
 
      check_validity();
    }

    point param_siz::get_recon_downsampling(ui32 comp_num) const
    {
      assert(comp_num < get_num_components());
 
      point factor(1u << skipped_resolutions, 1u << skipped_resolutions);
      const param_cod* cdp = cod->get_coc(comp_num);
      if (dfs && cdp && cdp->is_dfs_defined()) {
        ui16 dfs_idx = cdp->get_dfs_index();
        const param_dfs* d = dfs->get_dfs(dfs_idx);
        if (d != NULL)
          factor = d->get_res_downsamp(skipped_resolutions);
      }
      factor.x *= (ui32)cptr[comp_num].XRsiz;
      factor.y *= (ui32)cptr[comp_num].YRsiz;
      return factor;
    }

    point param_siz::get_recon_size(ui32 comp_num) const
    {
      assert(comp_num < get_num_components());
 
      point factor = get_recon_downsampling(comp_num);
      point r;
      r.x = ojph_div_ceil(Xsiz, factor.x) - ojph_div_ceil(XOsiz, factor.x);
      r.y = ojph_div_ceil(Ysiz, factor.y) - ojph_div_ceil(YOsiz, factor.y);
      return r;
    }
 

    //
    //
    //
    //
    //

    bool param_cap::write(outfile_base *file)
    {
      //marker size excluding header
      Lcap = 8;
 
      ui16 buf2;
      ui32 buf4;
      bool result = true;
 
      buf2 = JP2K_MARKER::CAP;
      buf2 = swap_bytes_if_le(buf2);
      result &= file->write(&buf2, sizeof(ui16)) == sizeof(ui16);
      buf2 = swap_bytes_if_le(Lcap);
      result &= file->write(&buf2, sizeof(ui16)) == sizeof(ui16);
      buf4 = swap_bytes_if_le(Pcap);
      result &= file->write(&buf4, sizeof(ui32)) == sizeof(ui32);
 
      buf2 = swap_bytes_if_le(Ccap[0]);
      result &= file->write(&buf2, sizeof(ui16)) == sizeof(ui16);
 
      return result;
    }

    void param_cap::read(infile_base *file)
    {
      if (file->read(&Lcap, 2) != 2)
        OJPH_ERROR(0x00050061, "error reading CAP marker");
      Lcap = swap_bytes_if_le(Lcap);
      if (file->read(&Pcap, 4) != 4)
        OJPH_ERROR(0x00050062, "error reading CAP marker");
      Pcap = swap_bytes_if_le(Pcap);
      ui32 count = population_count(Pcap);
      if (Pcap & 0xFFFDFFFF)
        OJPH_ERROR(0x00050063,
          "error Pcap in CAP has options that are not supported");
      if ((Pcap & 0x00020000) == 0)
        OJPH_ERROR(0x00050064,
          "error Pcap should have its 15th MSB set, Pcap^15. "
          " This is not a JPH file");
      for (ui32 i = 0; i < count; ++i)
        if (file->read(Ccap+i, 2) != 2)
          OJPH_ERROR(0x00050065, "error reading CAP marker");
      if (Lcap != 6 + 2 * count)
        OJPH_ERROR(0x00050066, "error in CAP marker length");
    }

    //
    //
    //
    //
    //

    bool param_cod::is_reversible() const
    {
      if (SPcod.wavelet_trans <= 1)
        return get_wavelet_kern() == local::param_cod::DWT_REV53;
      else {
        assert(atk != NULL);
        return atk->is_reversible();
      }
    }

    bool param_cod::write(outfile_base *file)
    {
      assert(type == COD_MAIN);
 
      //marker size excluding header
      Lcod = 12;
      Lcod = (ui16)(Lcod + (Scod & 1 ? 1 + SPcod.num_decomp : 0));
 
      ui8  buf1;
      ui16 buf2;
      bool result = true;
 
      buf2 = JP2K_MARKER::COD;
      buf2 = swap_bytes_if_le(buf2);
      result &= file->write(&buf2, sizeof(ui16)) == sizeof(ui16);
      buf2 = swap_bytes_if_le(Lcod);
      result &= file->write(&buf2, sizeof(ui16)) == sizeof(ui16);
      buf1 = Scod;
      result &= file->write(&buf1, sizeof(ui8)) == sizeof(ui8);
      buf1 = SGCod.prog_order;
      result &= file->write(&buf1, sizeof(ui8)) == sizeof(ui8);
      buf2 = swap_bytes_if_le(SGCod.num_layers);
      result &= file->write(&buf2, sizeof(ui16)) == sizeof(ui16);
      buf1 = SGCod.mc_trans;
      result &= file->write(&buf1, sizeof(ui8)) == sizeof(ui8);
      buf1 = SPcod.num_decomp;
      result &= file->write(&buf1, sizeof(ui8)) == sizeof(ui8);
      buf1 = SPcod.block_width;
      result &= file->write(&buf1, sizeof(ui8)) == sizeof(ui8);
      buf1 = SPcod.block_height;
      result &= file->write(&buf1, sizeof(ui8)) == sizeof(ui8);
      buf1 = SPcod.block_style;
      result &= file->write(&buf1, sizeof(ui8)) == sizeof(ui8);
      buf1 = SPcod.wavelet_trans;
      result &= file->write(&buf1, sizeof(ui8)) == sizeof(ui8);
      if (Scod & 1)
        for (int i = 0; i <= SPcod.num_decomp; ++i)
        {
          buf1 = SPcod.precinct_size[i];
          result &= file->write(&buf1, sizeof(ui8)) == sizeof(ui8);
        }
 
      return result;
    }

    bool param_cod::write_coc(outfile_base *file, ui32 num_comps)
    {
      assert(type == COD_MAIN);
      bool result = true;
      param_cod *p = this->next;
      while (p)
      {
        if (p->comp_idx < num_comps)
          result &= p->internal_write_coc(file, num_comps);
        p = p->next;
      }
      return result;
    }

    bool param_cod::internal_write_coc(outfile_base *file, ui32 num_comps)
    {
      assert(type == COC_MAIN);
 
      //marker size excluding header
      Lcod = num_comps < 257 ? 9 : 10;
      Lcod = (ui16)(Lcod + (Scod & 1 ? 1 + SPcod.num_decomp : 0));
 
      ui8  buf1;
      ui16 buf2;
      bool result = true;
 
      buf2 = JP2K_MARKER::COC;
      buf2 = swap_bytes_if_le(buf2);
      result &= file->write(&buf2, sizeof(ui16)) == sizeof(ui16);
      buf2 = swap_bytes_if_le(Lcod);
      result &= file->write(&buf2, sizeof(ui16)) == sizeof(ui16);
      if (num_comps < 257)
      {
        buf1 = (ui8)comp_idx;
        result &= file->write(&buf1, sizeof(ui8)) == sizeof(ui8);
      }
      else
      {
        buf2 = swap_bytes_if_le(comp_idx);
        result &= file->write(&buf2, sizeof(ui16)) == sizeof(ui16);
      }
      buf1 = Scod;
      result &= file->write(&buf1, sizeof(ui8)) == sizeof(ui8);
      buf1 = SPcod.num_decomp;
      result &= file->write(&buf1, sizeof(ui8)) == sizeof(ui8);
      buf1 = SPcod.block_width;
      result &= file->write(&buf1, sizeof(ui8)) == sizeof(ui8);
      buf1 = SPcod.block_height;
      result &= file->write(&buf1, sizeof(ui8)) == sizeof(ui8);
      buf1 = SPcod.block_style;
      result &= file->write(&buf1, sizeof(ui8)) == sizeof(ui8);
      buf1 = SPcod.wavelet_trans;
      result &= file->write(&buf1, sizeof(ui8)) == sizeof(ui8);
      if (Scod & 1)
        for (int i = 0; i <= SPcod.num_decomp; ++i)
        {
          buf1 = SPcod.precinct_size[i];
          result &= file->write(&buf1, sizeof(ui8)) == sizeof(ui8);
        }
 
      return result;
    }

    void param_cod::read(infile_base *file)
    {
      assert(type == COD_MAIN);
 
      if (file->read(&Lcod, 2) != 2)
        OJPH_ERROR(0x00050071, "error reading COD segment");
      Lcod = swap_bytes_if_le(Lcod);
      if (file->read(&Scod, 1) != 1)
        OJPH_ERROR(0x00050072, "error reading COD segment");
      if (file->read(&SGCod.prog_order, 1) != 1)
        OJPH_ERROR(0x00050073, "error reading COD segment");
      if (file->read(&SGCod.num_layers, 2) != 2)
      { OJPH_ERROR(0x00050074, "error reading COD segment"); }
      else
        SGCod.num_layers = swap_bytes_if_le(SGCod.num_layers);
      if (file->read(&SGCod.mc_trans, 1) != 1)
        OJPH_ERROR(0x00050075, "error reading COD segment");
      if (file->read(&SPcod.num_decomp, 1) != 1)
        OJPH_ERROR(0x00050076, "error reading COD segment");
      if (file->read(&SPcod.block_width, 1) != 1)
        OJPH_ERROR(0x00050077, "error reading COD segment");
      if (file->read(&SPcod.block_height, 1) != 1)
        OJPH_ERROR(0x00050078, "error reading COD segment");
      if (file->read(&SPcod.block_style, 1) != 1)
        OJPH_ERROR(0x00050079, "error reading COD segment");
      if (file->read(&SPcod.wavelet_trans, 1) != 1)
        OJPH_ERROR(0x0005007A, "error reading COD segment");
 
      if (get_num_decompositions() > 32
        || SPcod.block_width > 8
        || SPcod.block_height > 8
        || SPcod.block_width + SPcod.block_height > 8
        || (SPcod.block_style & 0x40) != 0x40
        || (SPcod.block_style & 0xB7) != 0x00)
        OJPH_ERROR(0x0005007D, "wrong settings in a COD-SPcod parameter");
      if ((SPcod.block_style & 0x40) != 0x40
        || (SPcod.block_style & 0xB7) != 0x00)
        OJPH_ERROR(0x0005007E, "unsupported settings in a COD-SPcod parameter");
 
      ui8 num_decompositions =  get_num_decompositions();
      if (Scod & 1) {
        for (int i = 0; i <= num_decompositions; ++i) {
          if (file->read(&SPcod.precinct_size[i], 1) != 1)
            OJPH_ERROR(0x0005007B, "error reading COD segment");
          if (i)
            if ((SPcod.precinct_size[i] & 0x0F) == 0 ||
              (SPcod.precinct_size[i] >> 4) == 0)
              OJPH_ERROR(0x0005007F,
                "Precinct width or height for resolutions other than the"
                " coarsest must be larger than 1; here, they are %d and %d,"
                " respectively.",
                1 << (SPcod.precinct_size[i] & 0x0F),
                1 << (SPcod.precinct_size[i] >> 4));
        }
      }
      if (Lcod != 12 + ((Scod & 1) ? 1 + SPcod.num_decomp : 0))
        OJPH_ERROR(0x0005007C, "error in COD segment length");
    }

    void param_cod::read_coc(infile_base* file, ui32 num_comps,
                             param_cod *top_cod)
    {
      assert(type == COC_MAIN);
      assert(top_cod != NULL);
 
      this->SGCod = top_cod->SGCod;
      this->top_cod = top_cod;
      if (file->read(&Lcod, 2) != 2)
        OJPH_ERROR(0x00050121, "error reading COC segment");
      Lcod = swap_bytes_if_le(Lcod);
      if (num_comps < 257) {
        ui8 t;
        if (file->read(&t, 1) != 1)
          OJPH_ERROR(0x00050122, "error reading COC segment");
        comp_idx = t;
      }
      else {
        if (file->read(&comp_idx, 2) != 2)
          OJPH_ERROR(0x00050123, "error reading COC segment");
        comp_idx = swap_bytes_if_le(comp_idx);
      }
      if (file->read(&Scod, 1) != 1)
        OJPH_ERROR(0x00050124, "error reading COC segment");
      if (Scod & 0xF8)
        OJPH_WARN(0x00050011,
          "Unsupported options in Scoc field of the COC segment");
      if (file->read(&SPcod.num_decomp, 1) != 1)
        OJPH_ERROR(0x00050125, "error reading COC segment");
      if (file->read(&SPcod.block_width, 1) != 1)
        OJPH_ERROR(0x00050126, "error reading COC segment");
      if (file->read(&SPcod.block_height, 1) != 1)
        OJPH_ERROR(0x00050127, "error reading COC segment");
      if (file->read(&SPcod.block_style, 1) != 1)
        OJPH_ERROR(0x00050128, "error reading COC segment");
      if (file->read(&SPcod.wavelet_trans, 1) != 1)
        OJPH_ERROR(0x00050129, "error reading COC segment");
 
      if (get_num_decompositions() > 32
        || SPcod.block_width > 8
        || SPcod.block_height > 8
        || SPcod.block_width + SPcod.block_height > 8
        || (SPcod.block_style & 0x40) != 0x40
        || (SPcod.block_style & 0xB7) != 0x00)
        OJPH_ERROR(0x0005012C, "wrong settings in a COC-SPcoc parameter");
      if ((SPcod.block_style & 0x40) != 0x40
        || (SPcod.block_style & 0xB7) != 0x00)
        OJPH_ERROR(0x0005012D, "unsupported settings in a COC-SPcoc parameter");
 
      ui8 num_decompositions =  get_num_decompositions();
      if (Scod & 1) {
        for (int i = 0; i <= num_decompositions; ++i) {
          if (file->read(&SPcod.precinct_size[i], 1) != 1)
            OJPH_ERROR(0x0005012A, "error reading COC segment");
          if (i)
            if ((SPcod.precinct_size[i] & 0x0F) == 0 ||
              (SPcod.precinct_size[i] >> 4) == 0)
              OJPH_ERROR(0x0005012E,
                "Precinct width or height for resolutions other than the"
                " coarsest must be larger than 1; here, they are %d and %d,"
                " respectively.",
                1 << (SPcod.precinct_size[i] & 0x0F),
                1 << (SPcod.precinct_size[i] >> 4));
        }
      }
      ui32 t = 9;
      t += num_comps < 257 ? 0 : 1;
      t += (Scod & 1) ? 1 + num_decompositions : 0;
      if (Lcod != t)
        OJPH_ERROR(0x0005012B, "error in COC segment length");
    }

    void param_cod::update_atk(param_atk* atk)
    {
      assert(type == COD_MAIN);
      this->atk = atk->get_atk(SPcod.wavelet_trans);
      if (this->atk == NULL)
        OJPH_ERROR(0x00050131, "A COD segment employs the DWT kernel "
          "atk = %d, but a corresponding ATK segment cannot be found.",
          SPcod.wavelet_trans);
      param_cod *p = next;
      while (p)
      {
        p->atk = atk->get_atk(p->SPcod.wavelet_trans);
        if (p->atk == NULL)
          OJPH_ERROR(0x00050132, "A COC segment employs the DWT kernel "
            "atk = %d, but a corresponding ATK segment cannot be found",
            SPcod.wavelet_trans);
        p = p->next;
      }
    }

    const param_cod* param_cod::get_coc(ui32 comp_idx) const
    {
      assert(this->type == COD_MAIN || this->top_cod->type == COD_MAIN);
      const param_cod *p, *q;
      if (this->type == COD_MAIN)
        q = p = this;
      else
        q = p = this->top_cod;
      while (p && p->comp_idx != comp_idx)
        p = p->next;
      return p ? p : q;
    }

    param_cod* param_cod::get_coc(ui32 comp_idx)
    {
      // cast object to constant
      const param_cod* const_p = const_cast<const param_cod*>(this);
      // call using the constant object, then cast to non-const
      return const_cast<param_cod*>(const_p->get_coc(comp_idx));
    }

    param_cod* param_cod::add_coc_object(ui32 comp_idx)
    {
      assert(type == COD_MAIN);
      param_cod *p = this;
      while (p->next != NULL)
        p = p->next;
      if (avail)
      {
        p->next = avail;
        avail = avail->next;
        p->next->init(this, (ui16)comp_idx);
      }
      else
        p->next = new param_cod(this, (ui16)comp_idx);
      return p->next;
    }

    param_cod* param_cod::get_or_add_coc(ui32 comp_idx)
    {
      assert(type == COD_MAIN);
      local::param_cod *p = get_coc(comp_idx);
      if (p == this)
        p = add_coc_object(comp_idx);
      return p;
    }

    //
    //
    //
    //
    //

    void param_qcd::check_validity(const param_siz& siz, const param_cod& cod)
    {
      assert(this->type == QCD_MAIN);
 
      ui32 num_comps = siz.get_num_components();
      trim_non_existing_components(num_comps);
 
      // initialize QCD based on the first component that is (a) associated with
      // COD and (b) does not have a COC, or the first component othewise.
      ui32 qcd_comp = 0;
      for (ui32 c = 0; c < num_comps; ++c)
      {
        if (cod.get_coc(c) == &cod && get_qcc(c) == this)
        {
          qcd_comp = c;
          break;
        }
      }
 
      // check if only the top QCD has qfactor set, if so, check if any
      // of the first component has COC and qfactor set properly
      if (this->qfactor != QFACTOR_UNSET)
      {
        if (num_comps < 3) // one or two components
        {
          for (ui32 i = 0; i < num_comps; ++i)
          {
            param_qcd* q = get_qcc(i);
            if (q == this)
            {
              q = add_qcc_object(i);
              set_qfactor(i, comp_type::OJPH_COMP_Y, this->qfactor);
            }
          }
        }
        else if (num_comps >= 3)
        {
          for (ui32 i = 0; i < num_comps; ++i) {
            param_qcd* q = get_qcc(i);
            if (q == this)
            {
              q = add_qcc_object(i);
              ui8 ci = (ui8)(i < 3u ? i : 0u);
              comp_type t = ojph::param_qcd::ui8_2_comp_type(ci);
              set_qfactor(i, t, this->qfactor);
            }
          }
        }
      }
 
      this->make_quant_steps(qcd_comp, cod, siz);
 
      // initialize every QCC, creating one for every component that (a) cannot
      // use QCD and (b) does not already have a QCC
      // NOTE: Qfactor always creates a QCC and QCD cannot be reused
      for (ui32 c = 0; c < num_comps; ++c)
      {
        param_qcd *qcc = this->get_qcc(c);
        const param_cod *coc = cod.get_coc(c);
 
        // check if a QCC exists for the component
        if (qcc == this)
        {
          // if none exists, do not create one if QCD can be reused
          if (!this->is_qcc_needed(c, *coc, siz))
            continue;
 
          qcc = this->add_qcc_object(c);
          qcc->set_delta(this->base_delta);
        }
 
        qcc->make_quant_steps(c, *coc, siz);
      }
    }

    void param_qcd::make_quant_steps(ui32 comp_num, const param_cod &cod,
                                     const param_siz &siz)
    {
      if (this->is_init)
        OJPH_ERROR(0x000501C1, "Quantization step sizes already initialized.");
 
      this->is_init = true;
 
      this->num_decomps = cod.get_num_decompositions();
      this->bit_depth = siz.get_bit_depth(comp_num);
      this->is_signed = siz.is_signed(comp_num);
      this->is_color_trans = cod.is_employing_color_transform();
      this->wavelet_kern = cod.get_wavelet_kern();
      this->sampling = siz.get_downsampling(comp_num);
      this->num_subbands = 1 + 3 * this->num_decomps;
 
      if (this->num_weights != 0)
      {
        if (this->wavelet_kern != param_cod::DWT_IRV97)
          OJPH_ERROR(0x000501C2, "Quantization weights for component %d can "
            "only be used with the irreversible transform.", comp_num);
        if (this->qfactor != QFACTOR_UNSET)
          OJPH_ERROR(0x000501C3, "Quantization weights for component %d "
            "cannot be used together with qfactor.", comp_num);
        if (this->num_weights != this->num_subbands)
          OJPH_ERROR(0x000501C4, "Component %d has %d quantization weights, "
            "but %d are needed for %d decomposition levels.", comp_num,
            this->num_weights, this->num_subbands, this->num_decomps);
      }
 
      if (this->wavelet_kern == param_cod::DWT_REV53)
        this->set_rev_quant(this->num_decomps, this->bit_depth,
                            comp_num < 3 ? this->is_color_trans : false);
      else if (this->wavelet_kern == param_cod::DWT_IRV97)
      {
        if (this->base_delta == -1.0f)
        {
          ui32 t = ojph_min(16, bit_depth);
          this->base_delta = 1.0f / (float)(1 << t);
        }
        else if (qfactor != QFACTOR_UNSET)
          OJPH_WARN(0x000501C5, "qstep for component %d is ignored, because "
            "qfactor is set.", comp_num);
 
        this->set_irrev_quant(this->num_decomps);
      }
    }

    bool param_qcd::is_qcc_needed(ui32 comp_num, const param_cod &cod,
                                  const param_siz &siz)
    {
      if (! this->is_init)
        OJPH_ERROR(0x000501D1, "Quantization step sizes not initialized.");
 
      return this->num_decomps != cod.get_num_decompositions() ||
              this->bit_depth != siz.get_bit_depth(comp_num) ||
              this->is_signed != siz.is_signed(comp_num) ||
              this->is_color_trans != cod.is_employing_color_transform() ||
              this->wavelet_kern != cod.get_wavelet_kern();
    }

    void param_qcd::set_qfactor(float qfactor) {
      assert(this->type == QCD_MAIN);
 
      if (qfactor < 1.0f || qfactor > 100.0f)
        OJPH_ERROR(0x00050181, "Qfactor must be between 1 and 100, "
          "but was set to %i.", qfactor);
 
      this->qfactor = qfactor;
    }

    void param_qcd::set_rev_quant(ui32 num_decomps, ui32 bit_depth,
                                  bool is_employing_color_transform)
    {
      ui32 B = bit_depth;
      B += is_employing_color_transform ? 1 : 0; //1 bit for RCT
      int s = 0;
      double bibo_l = bibo_gains::get_bibo_gain_l(num_decomps, true);
      ui32 X = (ui32) ceil(log(bibo_l * bibo_l) / M_LN2);
      SPqcd.u8[s++] = (ui8)(B + X);
      ui32 max_B_plus_X = (ui32)(B + X);
      for (ui32 d = num_decomps; d > 0; --d)
      {
        double bibo_l = bibo_gains::get_bibo_gain_l(d, true);
        double bibo_h = bibo_gains::get_bibo_gain_h(d - 1, true);
        X = (ui32) ceil(log(bibo_h * bibo_l) / M_LN2);
        SPqcd.u8[s++] = (ui8)(B + X);
        max_B_plus_X = ojph_max(max_B_plus_X, B + X);
        SPqcd.u8[s++] = (ui8)(B + X);
        max_B_plus_X = ojph_max(max_B_plus_X, B + X);
        X = (ui32) ceil(log(bibo_h * bibo_h) / M_LN2);
        SPqcd.u8[s++] = (ui8)(B + X);
        max_B_plus_X = ojph_max(max_B_plus_X, B + X);
      }
 
      if (max_B_plus_X > 38)
        OJPH_ERROR(0x00050151, "The specified combination of bit_depth, "
         "colour transform, and type of wavelet transform requires more than "
         "38 bits; it requires %d bits. This is beyond what is allowed in "
         "the JPEG2000 image coding format.", max_B_plus_X);
 
      int guard_bits = ojph_max(1, (si32)max_B_plus_X - 31);
      Sqcd = (ui8)(guard_bits << 5);
      s = 0;
      SPqcd.u8[s] = encode_SPqcd((ui8)(SPqcd.u8[s] - guard_bits));
      s++;
      for (ui32 d = num_decomps; d > 0; --d)
      {
        SPqcd.u8[s] = encode_SPqcd((ui8)(SPqcd.u8[s] - guard_bits));
        s++;
        SPqcd.u8[s] = encode_SPqcd((ui8)(SPqcd.u8[s] - guard_bits));
        s++;
        SPqcd.u8[s] = encode_SPqcd((ui8)(SPqcd.u8[s] - guard_bits));
        s++;
      }
    }

    void param_qcd::set_irrev_quant(ui32 num_decomps)
    {
      int guard_bits = 1;
      Sqcd = (ui8)((guard_bits<<5)|0x2); //one guard bit, scalar quantization
 
      float g_c = 1.0f;
      float delta_ref = base_delta;
      float power = 1.0f;
      const float* weights = visual_weights::get_no_weights();
 
      if (qfactor != QFACTOR_UNSET)
      {
        visual_weights::colour_format format =
          visual_weights::get_format(this->sampling);
        if (format == visual_weights::VW_COLOUR_FORMAT_ERROR)
          OJPH_ERROR(0x00050161, "Qfactor can only be used on components "
            "with 4:4:4, 4:2:2 or 4:2:0 sampling");
        if (this->ctype == comp_type::OJPH_COMP_Y &&
          this->sampling.x != 1 && this->sampling.y != 1)
          OJPH_ERROR(0x00050162, "Qfactor can only be used for a Y or "
            "luminance component when it is not downsampled.");
 
        // calculate component gain
        g_c = visual_weights::get_gain(this->ctype);
 
        // calculate delta_ref & power
        delta_ref = visual_weights::get_delta_ref(qfactor, bit_depth, power);
 
        // find visual weight
        weights = visual_weights::get_weights(format, this->ctype);
      }
 
      // returns sqrt(w_b), the square root of the weight, for a decomposition
      // level and subband (0:LL, 1:HL, 2:LH, 3:HH); coefficients are ordered
      // {LH1, HL1, HH1, ..., LLN}
      auto get_sqrt_weight = [&](ui32 d, ui32 sb) -> float
      {
        if (num_weights == 0)
          return visual_weights::get_weight(weights, d, sb);
        if (sb == 0)
          return this->weights[num_weights - 1];
        static const ui32 pos[4] = { 0, 1, 0, 2 };
        return this->weights[(d - 1) * 3 + pos[sb]];
      };
 
      // LL band
      ui32 b = 0;
      float w_b;
      float gain_l = sqrt_energy_gains::get_gain_l(num_decomps, false);
      w_b = get_sqrt_weight(num_decomps, b);
      w_b = std::pow(w_b, power);
      encode_SPqcd(b++, delta_ref / (gain_l * gain_l * g_c * w_b));
 
      // LL, HL, LH, HH, HL, LH, HH...
      for (ui32 d = num_decomps; d > 0; --d)
      {
        // compute square root of the enery gain factor W_g
        float gain_l = sqrt_energy_gains::get_gain_l(d, false);
        float gain_h = sqrt_energy_gains::get_gain_h(d - 1, false);
 
        w_b = get_sqrt_weight(d, 1);
        w_b = std::pow(w_b, power);
        encode_SPqcd(b++, delta_ref / (gain_h * gain_l * g_c * w_b));
        w_b = get_sqrt_weight(d, 2);
        w_b = std::pow(w_b, power);
        encode_SPqcd(b++, delta_ref / (gain_l * gain_h * g_c * w_b));
        w_b = get_sqrt_weight(d, 3);
        w_b = std::pow(w_b, power);
        encode_SPqcd(b++, delta_ref / (gain_h * gain_h * g_c * w_b));
      }
    }

    void param_qcd::encode_SPqcd(ui32 subband_index, float delta)
    {
      int exp = 0, mantissa;
      while (delta < 1.0f)
      { exp++; delta *= 2.0f; }
      mantissa = (int)round(delta * (float)(1<<11)) - (1<<11);
      // with rounding, there is a risk that the mantissa becomes
      // equal to 1<<11
      mantissa = mantissa < (1<<11) ? mantissa : 0x7FF;
      SPqcd.u16[subband_index] = (ui16)((exp << 11) | mantissa);
    }

    ui32 param_qcd::get_MAGB() const
    {
      ui32 B = 0;
 
      const param_qcd *p = this;
      while (p)
      {
        //this can be written better, but it is only executed once
        // this assumes a bi-directional wavelet (conventional DWT)
        ui32 num_decomps = (p->num_subbands - 1) / 3;
 
        int irrev = p->Sqcd & 0x1F;
        if (irrev == 0) //reversible
          for (ui32 i = 0; i < p->num_subbands; ++i) {
            ui32 t = p->decode_SPqcd(p->SPqcd.u8[i]);
            t += p->get_num_guard_bits() - 1u;
            B = ojph_max(B, t);
          }
        else if (irrev == 2) //scalar expounded
          for (ui32 i = 0; i < p->num_subbands; ++i)
          {
            ui32 nb = num_decomps - (i ? (i - 1) / 3 : 0); //decompsition level
            ui32 t = (p->SPqcd.u16[i] >> 11) + p->get_num_guard_bits() - nb;
            B = ojph_max(B, t);
          }
        else
          assert(0);
 
        p = p->next;
      }
 
      return B;
    }

    float param_qcd::get_irrev_delta(const param_dfs* dfs,
                                     ui32 num_decompositions, ui32 comp_num,
                                     ui32 resolution, ui32 subband) const
    {
      float arr[] = { 1.0f, 2.0f, 2.0f, 4.0f };
      if ((Sqcd & 0x1F) != 2)
        OJPH_ERROR(0x00050101, "There is something wrong in the configuration "
          "of the codestream; for component %d, the codestream defines an "
          "irreversible transform, for which the codestream provides a "
          "reversible (no quantization) step sizes in Sqcd/Sqcc.", comp_num);
 
      ui32 idx;
      if (dfs != NULL && dfs->exists())
        idx = dfs->get_subband_idx(num_decompositions, resolution, subband);
      else
        idx = resolution ? (resolution - 1) * 3 + subband : 0;
      if (idx >= num_subbands) {
        OJPH_INFO(0x00050102, "Trying to access quantization step size for "
          "subband %d when the QCD/QCC marker segment specifies "
          "quantization step sizes for %d subbands only.  To continue "
          "decoding, we are using the step size for subband %d, which can "
          "produce incorrect results",
          idx + 1, num_subbands, num_subbands - 1);
        idx = num_subbands - 1;
      }
      int eps = SPqcd.u16[idx] >> 11;
      float mantissa;
      mantissa = (float)((SPqcd.u16[idx] & 0x7FF) | 0x800) * arr[subband];
      mantissa /= (float)(1 << 11);
      mantissa /= (float)(1u << eps);
      return mantissa;
    }

    ui32 param_qcd::propose_precision(const param_cod* cod) const
    {
      ui32 comp_idx = cod->get_comp_idx();
      ui32 precision = 0;
      const param_cod *main =
        cod->get_coc(param_cod::OJPH_COD_DEFAULT);
      if (main->is_employing_color_transform() && comp_idx < 3)
      {
        for (ui32 i = 0; i < 3; ++i) {
          const param_qcd* p = this->get_qcc(i);
          precision = ojph_max(precision, p->get_largest_Kmax());
        }
      }
      else {
        precision = get_largest_Kmax();
      }
      // ``precision'' now holds the largest K_max, which excludes the sign
      // bit.
      // + 1 for the sign bit
      // + 1 because my block decoder/encoder does not supports up to 30
      //     bits (not 31), so we bump it by one more bit.
      return precision + 1 + 1;
    }

    ui32 param_qcd::get_num_guard_bits() const
    {
      return (Sqcd >> 5);
    }

    ui32 param_qcd::get_Kmax(const param_dfs* dfs, ui32 num_decompositions,
                             ui32 resolution, ui32 subband) const
    {
      ui32 idx;
      if (dfs != NULL && dfs->exists())
        idx = dfs->get_subband_idx(num_decompositions, resolution, subband);
      else
        idx = resolution ? (resolution - 1) * 3 + subband : 0;
      if (idx >= num_subbands) {
        OJPH_INFO(0x00050111, "Trying to access quantization step size for "
          "subband %d when the QCD/QCC marker segment specifies "
          "quantization step sizes for %d subbands only.  To continue "
          "decoding, we are using the step size for subband %d, which can "
          "produce incorrect results",
          idx + 1, num_subbands, num_subbands - 1);
        idx = num_subbands - 1;
      }
 
      int irrev = Sqcd & 0x1F;
      ui32 num_bits = 0;
      if (irrev == 0) // reversible; this is (10.22) from the J2K book
      {
        num_bits = decode_SPqcd(SPqcd.u8[idx]);
        num_bits = num_bits == 0 ? 0 : num_bits - 1;
      }
      else if (irrev == 1)
        assert(0);
      else if (irrev == 2) //scalar expounded
        num_bits = (SPqcd.u16[idx] >> 11) - 1;
      else
        assert(0);
 
      return num_bits + get_num_guard_bits();
    }

    ui32 param_qcd::get_largest_Kmax() const
    {
      int irrev = Sqcd & 0x1F;
      ui32 num_bits = 0;
      if (irrev == 0) // reversible; this is (10.22) from the J2K book
      {
        for (ui32 i = 0; i < num_subbands; ++i) {
          ui32 t = decode_SPqcd(SPqcd.u8[i]);
          num_bits = ojph_max(num_bits, t == 0 ? 0 : t - 1);
        }
      }
      else if (irrev == 1)
        assert(0);
      else if (irrev == 2) //scalar expounded
      {
        for (ui32 i = 0; i < num_subbands; ++i) {
          ui32 t = (SPqcd.u16[i] >> 11) - 1;
          num_bits = ojph_max(num_bits, t);
        }
      }
      else
        assert(0);
 
      return num_bits + get_num_guard_bits();
    }

    bool param_qcd::write(outfile_base *file)
    {
      int irrev = Sqcd & 0x1F;
 
      //marker size excluding header
      Lqcd = 3;
      if (irrev == 0)
        Lqcd = (ui16)(Lqcd + num_subbands);
      else if (irrev == 2)
        Lqcd = (ui16)(Lqcd + 2 * num_subbands);
      else
        assert(0);
 
      ui8  buf1;
      ui16 buf2;
      bool result = true;
 
      buf2 = JP2K_MARKER::QCD;
      buf2 = swap_bytes_if_le(buf2);
      result &= file->write(&buf2, sizeof(ui16)) == sizeof(ui16);
      buf2 = swap_bytes_if_le(Lqcd);
      result &= file->write(&buf2, sizeof(ui16)) == sizeof(ui16);
      buf1 = Sqcd;
      result &= file->write(&buf1, sizeof(ui8)) == sizeof(ui8);
 
      if (irrev == 0)
        for (ui32 i = 0; i < num_subbands; ++i)
        {
          buf1 = SPqcd.u8[i];
          result &= file->write(&buf1, sizeof(ui8)) == sizeof(ui8);
        }
      else if (irrev == 2)
        for (ui32 i = 0; i < num_subbands; ++i)
        {
          buf2 = swap_bytes_if_le(SPqcd.u16[i]);
          result &= file->write(&buf2, sizeof(ui16)) == sizeof(ui16);
        }
      else
        assert(0);
 
      return result;
    }

    bool param_qcd::write_qcc(outfile_base *file, ui32 num_comps)
    {
      assert(type == QCD_MAIN);
      bool result = true;
      param_qcd *p = this->next;
      while (p)
      {
        if (p->enabled)
          result &= p->internal_write_qcc(file, num_comps);
        p = p->next;
      }
      return result;
    }

    bool param_qcd::internal_write_qcc(outfile_base *file, ui32 num_comps)
    {
      int irrev = Sqcd & 0x1F;
 
      //marker size excluding header
      Lqcd = (ui16)(4 + (num_comps < 257 ? 0 : 1));
      if (irrev == 0)
        Lqcd = (ui16)(Lqcd + num_subbands);
      else if (irrev == 2)
        Lqcd = (ui16)(Lqcd + 2 * num_subbands);
      else
        assert(0);
 
      ui8  buf1;
      ui16 buf2;
      bool result = true;
 
      buf2 = JP2K_MARKER::QCC;
      buf2 = swap_bytes_if_le(buf2);
      result &= file->write(&buf2, sizeof(ui16)) == sizeof(ui16);
      buf2 = swap_bytes_if_le(Lqcd);
      result &= file->write(&buf2, sizeof(ui16)) == sizeof(ui16);
      if (num_comps < 257)
      {
        buf1 = (ui8)comp_idx;
        result &= file->write(&buf1, sizeof(ui8)) == sizeof(ui8);
      }
      else
      {
        buf2 = swap_bytes_if_le(comp_idx);
        result &= file->write(&buf2, sizeof(ui16)) == sizeof(ui16);
      }
      buf1 = Sqcd;
      result &= file->write(&buf1, sizeof(ui8)) == sizeof(ui8);
      if (irrev == 0)
        for (ui32 i = 0; i < num_subbands; ++i)
        {
          buf1 = SPqcd.u8[i];
          result &= file->write(&buf1, sizeof(ui8)) == sizeof(ui8);
        }
      else if (irrev == 2)
        for (ui32 i = 0; i < num_subbands; ++i)
        {
          buf2 = swap_bytes_if_le(SPqcd.u16[i]);
          result &= file->write(&buf2, sizeof(ui16)) == sizeof(ui16);
        }
      else
        assert(0);
 
      return result;
    }

    void param_qcd::trim_non_existing_components(ui32 num_comps)
    {
      assert(type == QCD_MAIN && comp_idx == OJPH_QCD_DEFAULT);
      param_qcd *p = this->next;
      while (p)
      {
        assert(p->type == QCC_MAIN);
        p->enabled = p->comp_idx < num_comps;
        p = p->next;
      }
    }

    void param_qcd::read(infile_base *file)
    {
      if (file->read(&Lqcd, 2) != 2)
        OJPH_ERROR(0x00050081, "error reading QCD marker");
      Lqcd = swap_bytes_if_le(Lqcd);
      if (file->read(&Sqcd, 1) != 1)
        OJPH_ERROR(0x00050082, "error reading QCD marker");
      if ((Sqcd & 0x1F) == 0)
      {
        num_subbands = (Lqcd - 3);
        if (num_subbands == 0)
          OJPH_ERROR(0x0005008A, "QCD marker segment that specifies no "
            "quantization informtion");
        if (num_subbands > MAX_SUBBAND_COUNT || Lqcd != 3 + num_subbands)
          OJPH_ERROR(0x00050083, "wrong Lqcd value of %d in QCD marker", Lqcd);
        for (ui32 i = 0; i < num_subbands; ++i)
          if (file->read(&SPqcd.u8[i], 1) != 1)
            OJPH_ERROR(0x00050084, "error reading QCD marker");
      }
      else if ((Sqcd & 0x1F) == 1)
      {
        num_subbands = 0;
        OJPH_ERROR(0x00050089,
          "Scalar derived quantization is not supported yet in QCD marker");
        if (Lqcd != 5)
          OJPH_ERROR(0x00050085, "wrong Lqcd value in QCD marker");
      }
      else if ((Sqcd & 0x1F) == 2)
      {
        num_subbands = (Lqcd - 3) / 2;
        if (num_subbands == 0)
          OJPH_ERROR(0x0005008B, "QCD marker segment that specifies no "
            "quantization informtion");
        if (num_subbands > MAX_SUBBAND_COUNT || Lqcd != 3 + 2 * num_subbands)
          OJPH_ERROR(0x00050086, "wrong Lqcd value of %d in QCD marker", Lqcd);
        for (ui32 i = 0; i < num_subbands; ++i)
        {
          if (file->read(&SPqcd.u16[i], 2) != 2)
            OJPH_ERROR(0x00050087, "error reading QCD marker");
          SPqcd.u16[i] = swap_bytes_if_le(SPqcd.u16[i]);
        }
      }
      else
        OJPH_ERROR(0x00050088, "wrong Sqcd value in QCD marker");
    }

    void param_qcd::read_qcc(infile_base *file, ui32 num_comps)
    {
      if (file->read(&Lqcd, 2) != 2)
        OJPH_ERROR(0x000500A1, "error reading QCC marker");
      Lqcd = swap_bytes_if_le(Lqcd);
      if (num_comps < 257)
      {
        ui8 v;
        if (file->read(&v, 1) != 1)
          OJPH_ERROR(0x000500A2, "error reading QCC marker");
        comp_idx = v;
      }
      else
      {
        if (file->read(&comp_idx, 2) != 2)
          OJPH_ERROR(0x000500A3, "error reading QCC marker");
        comp_idx = swap_bytes_if_le(comp_idx);
      }
      if (file->read(&Sqcd, 1) != 1)
        OJPH_ERROR(0x000500A4, "error reading QCC marker");
      ui32 offset = num_comps < 257 ? 4 : 5;
      if ((Sqcd & 0x1F) == 0)
      {
        num_subbands = (Lqcd - offset);
        if (num_subbands == 0)
          OJPH_ERROR(0x000500AC, "QCC marker segment that specifies no "
            "quantization informtion");
        if (num_subbands > MAX_SUBBAND_COUNT || Lqcd != offset + num_subbands)
          OJPH_ERROR(0x000500A5, "wrong Lqcd value of %d in QCC marker", Lqcd);
        for (ui32 i = 0; i < num_subbands; ++i)
          if (file->read(&SPqcd.u8[i], 1) != 1)
            OJPH_ERROR(0x000500A6, "error reading QCC marker");
      }
      else if ((Sqcd & 0x1F) == 1)
      {
        num_subbands = 0;
        OJPH_ERROR(0x000500AB,
          "Scalar derived quantization is not supported yet in QCC marker");
        if (Lqcd != offset)
          OJPH_ERROR(0x000500A7, "wrong Lqcc value in QCC marker");
      }
      else if ((Sqcd & 0x1F) == 2)
      {
        num_subbands = (Lqcd - offset) / 2;
        if (num_subbands == 0)
          OJPH_ERROR(0x000500AD, "QCC marker segment that specifies no "
            "quantization informtion");
        if (num_subbands>MAX_SUBBAND_COUNT || Lqcd != offset + 2*num_subbands)
          OJPH_ERROR(0x000500A8, "wrong Lqcc value of %d in QCC marker", Lqcd);
        for (ui32 i = 0; i < num_subbands; ++i)
        {
          if (file->read(&SPqcd.u16[i], 2) != 2)
            OJPH_ERROR(0x000500A9, "error reading QCC marker");
          SPqcd.u16[i] = swap_bytes_if_le(SPqcd.u16[i]);
        }
      }
      else
        OJPH_ERROR(0x000500AA, "wrong Sqcc value in QCC marker");
    }

    void param_qcd::set_delta(ui32 comp_idx, float delta)
    {
      assert(type == QCD_MAIN);
      param_qcd *p = get_qcc(comp_idx);
      if (p == this)
        p = add_qcc_object(comp_idx);
      p->set_delta(delta);
    }

    void param_qcd::set_weights(size_t num_weights, const float *weights)
    {
      if (weights == NULL || num_weights==0 || num_weights > MAX_SUBBAND_COUNT)
        OJPH_ERROR(0x000501A1, "The number of quantization weights must be "
          "between 1 and %d, but %zu were provided.",
          MAX_SUBBAND_COUNT, num_weights);
      for (size_t i = 0; i < num_weights; ++i)
        if (!(weights[i] > 0.0f))
          OJPH_ERROR(0x000501A2, "Quantization weights must be positive.");
      memcpy(this->weights, weights, num_weights * sizeof(float));
      this->num_weights = (ui32)num_weights;
    }

    void param_qcd::set_delta(ui32 comp_idx, float delta,
                              size_t num_weights, const float *weights)
    {
      assert(type == QCD_MAIN);
      param_qcd *p = get_qcc(comp_idx);
      if (p == this)
        p = add_qcc_object(comp_idx);
      p->set_delta(delta);
      p->set_weights(num_weights, weights);
    }

    void param_qcd::set_qfactor(ui32 comp_idx, comp_type ctype, float qfactor)
    {
      assert(this->type == QCD_MAIN);
 
      if (qfactor < 1.0f || qfactor > 100.0f)
        OJPH_ERROR(0x00050191, "Qfactor must be between 1 and 100, "
          "but was set to %i.", qfactor);
 
      param_qcd *p = get_qcc(comp_idx);
      if (p == this)
        p = add_qcc_object(comp_idx);
 
      p->qfactor = qfactor;
      p->ctype = ctype;
    }

    param_qcd* param_qcd::get_qcc(ui32 comp_idx)
    {
      // cast object to constant
      const param_qcd* const_p = const_cast<const param_qcd*>(this);
      // call using the constant object, then cast to non-const
      return const_cast<param_qcd*>(const_p->get_qcc(comp_idx));
    }

    const param_qcd* param_qcd::get_qcc(ui32 comp_idx) const
    {
      assert(this->type == QCD_MAIN || this->top_qcd->type == QCD_MAIN);
      const param_qcd *p, *q;
      if (this->type == QCD_MAIN)
        q = p = this;
      else
        q = p = this->top_qcd;
      while (p && p->comp_idx != comp_idx)
        p = p->next;
      return p ? p : q;
    }

    param_qcd* param_qcd::add_qcc_object(ui32 comp_idx)
    {
      assert(type == QCD_MAIN);
      param_qcd *p = this;
      while (p->next != NULL)
        p = p->next;
      if (avail)
      {
        p->next = avail;
        avail = avail->next;
        p->next->init(this, (ui16)comp_idx);
      }
      else
        p->next = new param_qcd(this, (ui16)comp_idx);
      return p->next;
    }

    //
    //
    //
    //
    //

    void nlt_rec::prepare_for_decoding()
    {
      double d = 1.0 / (double)((1ull << 32) - 1);
      fd_min = (float)((double)d_min * d);
      fd_max = (float)((double)d_max * d);
      delta = (fd_max - fd_min) / (float)(num_points - 1);
      inv_delta = (float)(num_points - 1) / (fd_max - fd_min);
      multiplier = (float)(1ull << get_bit_depth());
      float divider = (float)(1.0 / (double)(1ull << get_bit_depth()));
      if (bytes_per_point == 1) {
        ui8* sp = (ui8*)marker_points;
        float* dp = dec_points;
        for (ui32 i = 0; i < num_points; ++i)
          *dp++ = *sp++ * divider;
        dec_points[-1] = dec_points[0];
        dec_points[num_points] = dec_points[num_points - 1];
      }
      else if (bytes_per_point == 2) {
        ui16* sp = (ui16*)marker_points;
        float* dp = dec_points;
        for (ui32 i = 0; i < num_points; ++i)
          *dp++ = *sp++ * divider;
        dec_points[-1] = dec_points[0];
        dec_points[num_points] = dec_points[num_points - 1];
      }
      else if (bytes_per_point == 4) {
        ui32* sp = (ui32*)marker_points;
        float* dp = dec_points;
        for (ui32 i = 0; i < num_points; ++i)
          *dp++ = (float)*sp++ * divider;
        dec_points[-1] = dec_points[0];
        dec_points[num_points] = dec_points[num_points - 1];
      }
      else
        assert(0);
    }

    template<typename T>
    static inline void
    prepare_for_encoding_approx_for_type(const void* marker_points, float div,
      float& ft_min, float& ft_max, float& delta, float& inv_delta,
      float fd_min, float fd_max, ui32 num_points, ui32 enc_num_points,
      float* approx_enc_points)
    {
      const T* p = (const T*)marker_points;
      ft_min = (float)p[0] * div;
      ft_max = (float)p[num_points - 1] * div;
      delta = (ft_max - ft_min) / (float)(enc_num_points - 1);
      inv_delta = (float)(enc_num_points - 1) / (ft_max - ft_min);
 
      ui32 k = 0;
      float y_k = (float)p[k] * div, y_kp1 = (float)p[k + 1] * div;
      float dt = (fd_max - fd_min) / (float)(num_points - 1);
      float d_k = fd_min, d_kp1 = fd_min + dt;
      for (ui32 i = 1; i < enc_num_points - 1; ++i)
      {
        float z = ft_min + (float)i * delta;
        while (k + 1 < num_points - 1 && z >= y_kp1)
        {
          ++k;
          d_k   = d_kp1;
          d_kp1 = fd_min + (float)(k + 1) * dt;
          y_k   = y_kp1;
          y_kp1 = (float)p[k + 1] * div;
        }
        approx_enc_points[i] = d_k + (z - y_k) * dt / (y_kp1 - y_k);
      }
    }

    void nlt_rec::prepare_for_encoding_approx()
    {
      double d = 1.0 / (double)((1ull << 32) - 1);
      fd_min = (float)((double)d_min * d);
      fd_max = (float)((double)d_max * d);
 
      // enc_points holds positions in the domain of the LUT, whose first
      // and last entries are fd_min and fd_max, not LUT values
      approx_enc_points[-2] = fd_min;
      approx_enc_points[-1] = fd_min;
      approx_enc_points[ 0] = fd_min;
      approx_enc_points[enc_num_points - 1] = fd_max;
      approx_enc_points[enc_num_points] = fd_max;
      approx_enc_points[enc_num_points + 1] = fd_max;
 
      // create lookup table for encoding
      float mul = (float)((1ull << pt_val) - 1);
      float div = 1.0f / mul;
      if (bytes_per_point == 1)
        prepare_for_encoding_approx_for_type<ui8>(marker_points, div,
          ft_min, ft_max, delta, inv_delta,
          fd_min, fd_max, num_points, enc_num_points,
          approx_enc_points);
      else if (bytes_per_point == 2)
        prepare_for_encoding_approx_for_type<ui16>(marker_points, div,
          ft_min, ft_max, delta, inv_delta,
          fd_min, fd_max, num_points, enc_num_points,
          approx_enc_points);
      else if (bytes_per_point == 4)
        prepare_for_encoding_approx_for_type<ui32>(marker_points, div,
          ft_min, ft_max, delta, inv_delta,
          fd_min, fd_max, num_points, enc_num_points,
          approx_enc_points);
 
      multiplier = (float)(1.0 / (double)(1ull << get_bit_depth()));
    }

    void nlt_rec::prepare_for_encoding_precise()
    {
      double d = 1.0 / (double)((1ull << 32) - 1);
      fd_min = (float)((double)d_min * d);
      fd_max = (float)((double)d_max * d);
 
      float mul = (float)((1ull << pt_val) - 1);
      float div = 1.0f / mul;
 
      // convert decoder lookup table to floats
      if (bytes_per_point == 1) {
        const ui8* p = (const ui8*)marker_points;
        for (ui32 i = 0; i < num_points; ++i)
          precise_dec_lut[i] = (float)p[i] * div;
        ft_min = precise_dec_lut[0];
        ft_max = precise_dec_lut[num_points - 1];
      }
      else if (bytes_per_point == 2) {
        const ui16* p = (const ui16*)marker_points;
        for (ui32 i = 0; i < num_points; ++i)
          precise_dec_lut[i] = (float)p[i] * div;
        ft_min = precise_dec_lut[0];
        ft_max = precise_dec_lut[num_points - 1];
      }
      else if (bytes_per_point == 4) {
        const ui32* p = (const ui32*)marker_points;
        for (ui32 i = 0; i < num_points; ++i)
          precise_dec_lut[i] = (float)p[i] * div;
        ft_min = precise_dec_lut[0];
        ft_max = precise_dec_lut[num_points - 1];
      }
      precise_dec_lut[-2] = precise_dec_lut[0];
      precise_dec_lut[-1] = precise_dec_lut[0];
      precise_dec_lut[num_points]     = precise_dec_lut[num_points - 1];
      precise_dec_lut[num_points + 1] = precise_dec_lut[num_points - 1];
      delta = (fd_max - fd_min) / (float)(num_points - 1);
      inv_delta = (float)(enc_num_points - 1) / (ft_max - ft_min);
 
      multiplier = (float)(1.0 / (double)(1ull << get_bit_depth()));
 
      // find index into decoder lookup table
      float td = (ft_max - ft_min) / (float)(enc_num_points - 1);
      precise_max_steps = 0;
      ui16 k = 0, old_k = 0;
      for (ui32 i = 0; i < enc_num_points; ++i)
      {
        float z = ft_min + (float)i * td;
        while (k < num_points - 1 && z >= precise_dec_lut[k + 1])
          ++k;
        precise_enc_dec_indices[i] = k;
        precise_max_steps = ojph_max(precise_max_steps, (ui32)(k - old_k));
        old_k = k;
      }
    }

    void nlt_rec::prepare_for_encoding()
    {
      if (approx_encoding_nlt)
        prepare_for_encoding_approx();
      else
        prepare_for_encoding_precise();
    }

    int param_nlt::find_unsupported_nlt(const param_siz& siz,
                                        const param_cod& cod) const
    {
      ui32 num_comps = siz.get_num_components();
      for (ui32 c = 0; c < num_comps; ++c)
      {
        const nlt_rec* rec = get_nlt_rec(c);
        if (rec == NULL)
          continue;
        ui8 type = rec->get_type();
        // the LUT style nonlinearities are implemented for the irreversible
        // (9/7) wavelet only
        if ((type == nonlinearity::OJPH_NLT_LUT_STYLE_NLT ||
             type == nonlinearity::OJPH_NLT_BINARY_COMPLEMENT_PLUS_LUT) &&
            cod.get_coc(c)->is_reversible())
          return (int)c;
      }
      return -1;
    }

    void param_nlt::check_validity(param_siz& siz, const param_cod& cod)
    {
      if (is_any_enabled() == false)
        return;
 
      // the reversible wavelet cannot be used with a LUT style nonlinearity,
      // because applying that nonlinearity would need the samples of the
      // component to be transformed in a way this library cannot invert;
      // refusing it here is better than writing a codestream whose nonlinearity
      // a decoder cannot undo
      int comp = find_unsupported_nlt(siz, cod);
      if (comp >= 0)
        OJPH_ERROR(0x000501B1, "The LUT style nonlinearities (type 2, LUT "
          "style, and type 4, binary complement followed by a LUT) are "
          "supported with the irreversible (9/7) wavelet only; component %d of "
          "this codestream is coded with the reversible (5/3) wavelet. Use the "
          "irreversible wavelet for that component, use the binary complement "
          "nonlinearity (type 3), or do not use a nonlinearity with it.",
          comp);
 
      if (this->enabled && this->rec.Tnlt == nonlinearity::OJPH_NLT_UNDEFINED)
        this->enabled = false;
 
      if (this->enabled &&
          this->rec.Tnlt == nonlinearity::OJPH_NLT_BINARY_COMPLEMENT_NLT)
      {
        bool all_same = true;
        ui32 num_comps = siz.get_num_components();
 
        // first stage; find out if all components captured by the default
        // entry (ALL_COMPS) has the same bit_depth/signedness,
        // while doing this, set the BDnlt for components not captured by the
        // default entry (ALL_COMPS)
        ui32 bit_depth = 0;      // unknown yet
        bool is_signed = false;  // unknown yet
        for (ui32 c = 0; c < num_comps; ++c)
        { // captured by ALL_COMPS
          param_nlt* p = get_nlt_object(c);
          if (p == NULL || !p->enabled)
          {
            if (bit_depth != 0)
            {
              // we have seen an undefined component previously
              all_same = all_same && (bit_depth == siz.get_bit_depth(c));
              all_same = all_same && (is_signed == siz.is_signed(c));
            }
            else
            {
              // this is the first component which has not type 3 nlt definition
              bit_depth = siz.get_bit_depth(c);
              is_signed = siz.is_signed(c);
            }
          }
          else
          { // can be type 0 or type 3
            p->rec.BDnlt = (ui8)(siz.get_bit_depth(c) - 1);
            p->rec.BDnlt = (ui8)(p->rec.BDnlt | (siz.is_signed(c) ? 0x80 : 0));
          }
        }
 
        if (all_same && bit_depth != 0)
        { // all the same, and some components are captured by ALL_COMPS
          this->rec.BDnlt = (ui8)(bit_depth - 1);
          this->rec.BDnlt = (ui8)(this->rec.BDnlt | (is_signed ? 0x80 : 0));
        }
        else if (!all_same)
        { // have different settings or no component is captured by ALL_COMPS
          this->enabled = false;
          for (ui32 c = 0; c < num_comps; ++c)
          {
            param_nlt* p = get_nlt_object(c);
            if (p == NULL || !p->enabled)
            { // captured by ALL_COMPS
              if (p == NULL)
                p = add_object(c);
              p->enabled = true;
              p->rec.Tnlt = nonlinearity::OJPH_NLT_BINARY_COMPLEMENT_NLT;
              p->rec.BDnlt = (ui8)(siz.get_bit_depth(c) - 1);
              p->rec.BDnlt = (ui8)(p->rec.BDnlt | (siz.is_signed(c) ? 0x80:0));
            }
          }
        }
      }
      else {
        // fill NLT segment markers with correct information
        ui32 num_comps = siz.get_num_components();
        for (ui32 c = 0; c < num_comps; ++c)
        {
          param_nlt* p = get_nlt_object(c);
          if (p != NULL && p->enabled)
          { // can be type 0 or type 3
            p->rec.BDnlt = (ui8)(siz.get_bit_depth(c) - 1);
            p->rec.BDnlt = (ui8)(p->rec.BDnlt | (siz.is_signed(c) ? 0x80 : 0));
          }
        }
      }
 
      trim_non_existing_components(siz.get_num_components());
 
      if (is_any_enabled() == true)
        siz.set_Rsiz_flag(param_siz::RSIZ_EXT_FLAG | param_siz::RSIZ_NLT_FLAG);
    }

    void param_nlt::set_nonlinear_transform(ui32 comp_num, ui8 nl_type)
    {
      if (nl_type != ojph::param_nlt::OJPH_NLT_NO_NLT &&
          nl_type != ojph::param_nlt::OJPH_NLT_BINARY_COMPLEMENT_NLT)
      OJPH_ERROR(0x00050171, "Nonliearities other than type 0 "
        "(No Nonlinearity) or type 3 (Binary Binary Complement to Sign "
        "Magnitude Conversion) are not supported in this function call");
      param_nlt* p = get_nlt_object(comp_num);
      if (p == NULL)
        p = add_object(comp_num);
      p->Lnlt = 6;
      p->rec.Tnlt = nl_type;
      p->enabled = true;
    }

    void param_nlt::set_nonlinear_transform(ui32 comp_num,
                                            ui8 decoded_bit_depth,
                                            bool decoded_signedness,
                                            ui32 d_min, ui32 d_max, ui8 pt_val,
                                            ui16 num_points, void* points,
                                            ui8 nl_type,
                                            bool approx_encoding_nlt,
                                            ui32 num_enc_points)
    {
      if (nl_type != ojph::param_nlt::OJPH_NLT_LUT_STYLE_NLT &&
          nl_type != ojph::param_nlt::OJPH_NLT_BINARY_COMPLEMENT_PLUS_LUT)
        OJPH_ERROR(0x000501A1, "Nonliearities other than type 1 "
          "(LUT Nonlinearity) or type 4 (Binary Binary Complement to Sign "
          "Magnitude Conversion followed by LUT Nonlinearity) are not "
          "supported in this function call");
      if (pt_val == 0 || pt_val > 32)
        OJPH_ERROR(0x000501A2, "pt_val must be larger than 0 and "
          "smaller than or equal to 32; you provided, %d", pt_val);
      if (decoded_bit_depth == 0 || decoded_bit_depth > 38)
        OJPH_ERROR(0x000501A3, "decoded_bit_depth must be larger than 0 and "
          "smaller than or equal to 38; you provided, %d", decoded_bit_depth);
      if (d_min >= d_max)
        OJPH_ERROR(0x000501A4, "d_max must be larger than d_min; you "
          "provided %d and %d", d_min, d_max);
      if (num_points < 2 || num_points > 8192)
        OJPH_ERROR(0x000501A5, "The number of points must be larger or equal "
          "to 2 and smaller or equal to 8192 -- the code will convert this "
          "number to a number in the range 1 to 8191 for the NLT marker "
          "segement; you provided %d", num_points);
      if (num_enc_points == 1)
        OJPH_ERROR(0x000501AA, "The num_enc_points cannot be 1");
 
      param_nlt* p = get_nlt_object(comp_num);
      if (p == NULL)
        p = add_object(comp_num);
 
      p->Lnlt = (ui16)(17u + (ui32)num_points * (ui32)p->rec.get_bpp(pt_val));
      p->rec.BDnlt =
        (ui8)((decoded_bit_depth - 1) | (decoded_signedness ? 0x80 : 0));
      p->rec.Tnlt = nl_type;
 
      p->rec.d_min = d_min;
      p->rec.d_max = d_max;
      p->rec.pt_val = pt_val;
      p->rec.num_points = num_points;
      p->rec.bytes_per_point = p->rec.get_bpp(pt_val);
      p->rec.approx_encoding_nlt = approx_encoding_nlt;
 
      // Check that the LUT has increasing entries or has almost flat segments
      ui32 v_min = 0, v_max = 0;
      ui32 smallest_gap = UINT_MAX;
      if (p->rec.bytes_per_point == 1)
      {
        ui8* p = (ui8*)points;
        for (int k = 0; k < num_points - 1; ++k)
        {
          if (p[k+1] > p[k])
            smallest_gap = ojph_min(smallest_gap, (ui32)(p[k+1] - p[k]));
          else
            OJPH_ERROR(0x000501A6, "The LUT must have increasing sequence "
              "of numbers; here we have at LUT[%d+1] <= LUT[%d], with values "
              "%d and %d.", k, k, p[k+1], p[k]);
        }
        v_min = (ui32)p[0]; v_max = (ui32)p[num_points - 1];
      }
      else if (p->rec.bytes_per_point == 2)
      {
        ui16* p = (ui16*)points;
        for (int k = 0; k < num_points - 1; ++k)
        {
          if (p[k+1] > p[k])
            smallest_gap = ojph_min(smallest_gap, (ui32)(p[k+1] - p[k]));
          else
            OJPH_ERROR(0x000501A7, "The LUT must have increasing sequence "
              "of numbers; here we have at LUT[%d+1] <= LUT[%d], with values "
              "%d and %d.", k, k, p[k+1], p[k]);
        }
        v_min = (ui32)p[0]; v_max = (ui32)p[num_points - 1];
      }
      else if (p->rec.bytes_per_point == 4)
      {
        ui32* p = (ui32*)points;
        for (int k = 0; k < num_points - 1; ++k)
        {
          if (p[k+1] > p[k])
            smallest_gap = ojph_min(smallest_gap, p[k+1] - p[k]);
          else
            OJPH_ERROR(0x000501A8, "The LUT must have increasing sequence "
              "of numbers; here we have at LUT[%d+1] <= LUT[%d], with values "
              "%d and %d.", k, k, p[k+1], p[k]);
        }
        v_min = (ui32)p[0]; v_max = (ui32)p[num_points - 1];
      }
      else
        assert(0);
 
      // find ceil of the ratio to a power of 2
      ui32 ienc_pnts = num_enc_points;
      if (num_enc_points == 0)
      {
        float enc_pnts =
          std::ceil((float)(v_max-v_min) / (float)(smallest_gap));
        if (enc_pnts > 8192.0f)
        {
          ienc_pnts = 8192;
          OJPH_WARN(0x000501A1, "Encoding with LUT is performed with an "
            "encoding LUT, derived from the LUT you provided; however, "
            "because the provided LUT has almost flat segment or segments, "
            "these are hard to invert.  We are limiting the encoding "
            "LUT to 8192 entries, which means that some segment of the "
            "LUT table might be ignored during encoding.")
        }
        else {
          ienc_pnts = (ui32)enc_pnts;
          ienc_pnts = 32 - count_leading_zeros(ienc_pnts);
          ienc_pnts = 1u << ienc_pnts;
          ienc_pnts = ojph_max(ienc_pnts, 2);
        }
      }
 
      ui32 len = p->rec.cal_store_size_for_encoding(ienc_pnts);
      if (len > p->rec.store_size)
      {
        if (p->rec.points_store)
          delete[] (ui8*)p->rec.points_store;
        p->rec.store_size = 0;
        p->rec.points_store = new (std::nothrow) ui8[len];
        if (p->rec.points_store == NULL)
          OJPH_ERROR(0x000501A9, "Failed to allocated memory");
        p->rec.store_size = len;
      }
      p->rec.assign_pointers_for_encoding();
      memcpy(p->rec.marker_points, points, num_points * p->rec.bytes_per_point);
      p->rec.prepare_for_encoding();
 
      p->enabled = true;
    }

    bool
    param_nlt::get_nonlinear_transform(ui32 comp_num, ui8& bit_depth,
                                       bool& is_signed, ui8& nl_type) const
    {
      // top level must be ALL_COMPS; cannot call this functions from
      // non-top nlt
      assert(Cnlt == special_comp_num::ALL_COMPS);
      const param_nlt* p = get_nlt_object(comp_num);
      p = (p && p->enabled) ? p : this;
      if (p->enabled)
      {
        bit_depth = (ui8)((p->rec.BDnlt & 0x7F) + 1);
        bit_depth = bit_depth <= 38 ? bit_depth : 38;
        is_signed = (p->rec.BDnlt & 0x80) == 0x80;
        nl_type = (nonlinearity)p->rec.Tnlt;
        return true;
      }
      return false;
    }

    bool param_nlt::write(outfile_base* file) const
    {
      if (is_any_enabled() == false)
        return true;
 
      ui8 buf1;
      ui16 buf2;
      ui32 buf4;
      bool result = true;
      const param_nlt* p = this;
      while (p)
      {
        if (p->enabled)
        {
          buf2 = JP2K_MARKER::NLT;
          buf2 = swap_bytes_if_le(buf2);
          result &= file->write(&buf2, sizeof(ui16)) == sizeof(ui16);
          buf2 = swap_bytes_if_le(p->Lnlt);
          result &= file->write(&buf2, sizeof(ui16)) == sizeof(ui16);
          buf2 = swap_bytes_if_le(p->Cnlt);
          result &= file->write(&buf2, sizeof(ui16)) == sizeof(ui16);
          result &= file->write(&p->rec.BDnlt, 1) == 1;
          result &= file->write(&p->rec.Tnlt, 1) == 1;
          if (p->rec.Tnlt == nonlinearity::OJPH_NLT_LUT_STYLE_NLT ||
            p->rec.Tnlt == nonlinearity::OJPH_NLT_BINARY_COMPLEMENT_PLUS_LUT)
          {
            buf2 = (ui16)(p->rec.num_points - 1);
            buf2 = swap_bytes_if_le(buf2);
            result &= file->write(&buf2, sizeof(ui16)) == sizeof(ui16);
            buf4 = swap_bytes_if_le(p->rec.d_min);
            result &= file->write(&buf4, sizeof(ui32)) == sizeof(ui32);
            buf4 = swap_bytes_if_le(p->rec.d_max);
            result &= file->write(&buf4, sizeof(ui32)) == sizeof(ui32);
            // pt_val is a single byte in the marker segment; it must not be
            // written through the address of the 32 bit member, which would
            // write its most significant (and usually zero) byte on a
            // big-endian machine
            buf1 = (ui8)p->rec.pt_val;
            result &= file->write(&buf1, sizeof(ui8)) == sizeof(ui8);
 
            ui32 len = p->rec.bytes_per_point * p->rec.num_points;
            if (p->rec.bytes_per_point == 1)
              result &= file->write(p->rec.marker_points, len) == len;
            else if (p->rec.bytes_per_point == 2)
            {
              char* sp = (char*)p->rec.marker_points;
              for (ui32 i = 0; i < p->rec.num_points; ++i, sp += 2)
              {
                memcpy(&buf2, sp, 2); // not to violate strict aliasing
                buf2 = swap_bytes_if_le(buf2);
                result &= file->write(&buf2, sizeof(ui16)) == sizeof(ui16);
              }
            }
            else
            {
              char* sp = (char*)p->rec.marker_points;
              for (ui32 i = 0; i < p->rec.num_points; ++i, sp += 4)
              {
                memcpy(&buf4, sp, 4); // not to violate strict aliasing
                buf4 = swap_bytes_if_le(buf4);
                result &= file->write(&buf4, sizeof(ui32)) == sizeof(ui32);
              }
            }
          }
        }
        p = p->next;
      }
      return result;
    }

    bool param_nlt::read(infile_base* file)
    {
      ui16 buf2_len;
      ui16 buf2_comp;
      ui8  buf1_BDnlt;
      ui8  buf1_Tnlt;
      bool result = true;
 
      result &= file->read(&buf2_len, sizeof(ui16)) == sizeof(ui16);
      result &= file->read(&buf2_comp, sizeof(ui16)) == sizeof(ui16);
      result &= file->read(&buf1_BDnlt, sizeof(ui8)) == sizeof(ui8);
      result &= file->read(&buf1_Tnlt, sizeof(ui8)) == sizeof(ui8);
      if (result == false)
        OJPH_ERROR(0x000501E1, "Error reading NLT marker, Location 1");
 
      if (buf1_Tnlt != nonlinearity::OJPH_NLT_NO_NLT &&
        buf1_Tnlt != nonlinearity::OJPH_NLT_BINARY_COMPLEMENT_NLT &&
        buf1_Tnlt != nonlinearity::OJPH_NLT_LUT_STYLE_NLT &&
        buf1_Tnlt != nonlinearity::OJPH_NLT_BINARY_COMPLEMENT_PLUS_LUT)
        OJPH_ERROR(0x00050141, "Unsupported nonlinearity type in an "
          "NLT marker segment, with a value of", buf1_Tnlt);
 
      ui16 length = swap_bytes_if_le(buf2_len);
      if (buf1_Tnlt == nonlinearity::OJPH_NLT_NO_NLT ||
        buf1_Tnlt == nonlinearity::OJPH_NLT_BINARY_COMPLEMENT_NLT)
      {
        if (length != 6)
          OJPH_ERROR(0x00050142, "Somethins is wrong with this NLT "
            "marker segment -- it has an incorrect length of %d",
            length);
      }
 
      ui16 comp = swap_bytes_if_le(buf2_comp);
      if (comp > 16383 && comp != special_comp_num::ALL_COMPS)
        OJPH_ERROR(0x00050143, "Something is wrong with this NLT "
          "marker segment -- it has an incorrect component index of %d",
          comp);
 
      param_nlt* p = get_nlt_object(comp);
      if (p == NULL)
        p = add_object(comp);
      p->enabled = true;
      p->Cnlt = comp;
      p->rec.BDnlt = buf1_BDnlt;
      if ((buf1_BDnlt & 0x7F) + 1 > 38)
        OJPH_ERROR(0x00050144, "The NLT marker segment has an invalid bit "
          "depth of %d", (buf1_BDnlt & 0x7F) + 1);
      p->rec.Tnlt = buf1_Tnlt;
 
      if (buf1_Tnlt == nonlinearity::OJPH_NLT_LUT_STYLE_NLT ||
        buf1_Tnlt == nonlinearity::OJPH_NLT_BINARY_COMPLEMENT_PLUS_LUT)
      {
        ui16 buf2_num_points;
        ui32 buf4_d_min, buf4_d_max;
        ui8 buf_pt_val;
 
        result &= file->read(&buf2_num_points, sizeof(ui16)) == sizeof(ui16);
        result &= file->read(&buf4_d_min, sizeof(ui32)) == sizeof(ui32);
        result &= file->read(&buf4_d_max, sizeof(ui32)) == sizeof(ui32);
        result &= file->read(&buf_pt_val, sizeof(ui8)) == sizeof(ui8);
        if (result == false)
          OJPH_ERROR(0x000501E2, "Error reading NLT marker, Location 2");
 
        buf2_num_points = swap_bytes_if_le(buf2_num_points);
        buf4_d_min = swap_bytes_if_le(buf4_d_min);
        buf4_d_max = swap_bytes_if_le(buf4_d_max);
 
        if (buf2_num_points < 1 || buf2_num_points > 8191)
          OJPH_ERROR(0x00050145, "The NLT marker segment has an invalid "
            "number of points field", buf2_num_points);
        if (buf4_d_min >= buf4_d_max)
          OJPH_ERROR(0x00050146, "In an NLT marker segment, Dmin must smaller "
            "than Dmax; here we have %d and %d", buf4_d_min, buf4_d_max);
        if (buf_pt_val == 0 || buf_pt_val > 32)
          OJPH_ERROR(0x00050147, "In an NLT marker segment, the value of "
            "pt_val must be larger than 0 and smaller or equal to 32; "
            "here, we have %d", buf_pt_val);
 
        p->rec.num_points = buf2_num_points + 1;
        p->rec.d_min = buf4_d_min;
        p->rec.d_max = buf4_d_max;
        p->rec.pt_val = buf_pt_val;
        p->rec.bytes_per_point = p->rec.get_bpp(buf_pt_val);
 
        // read the points
        ui32 len = p->rec.cal_store_size_for_decoding();
        if (p->rec.store_size < len)
        {
          if (p->rec.points_store)
            delete[] (ui8*)p->rec.points_store;
          p->rec.store_size = 0;
          p->rec.points_store = new (std::nothrow) ui8[len];
          if (p->rec.points_store == NULL)
            OJPH_ERROR(0x00050148, "Failed to allocated memory");
          p->rec.store_size = len;
        }
        p->rec.assign_pointers_for_decoding();
 
        if (p->rec.bytes_per_point == 1)
          result &=
            file->read(p->rec.marker_points, p->rec.num_points)
              == p->rec.num_points;
        else if (p->rec.bytes_per_point == 2)
        {
          ui16 buf2;
          char* dp = (char*)p->rec.marker_points;
          for (ui32 i = 0; i < p->rec.num_points; ++i, dp += 2)
          {
            result &= file->read(&buf2, sizeof(ui16)) == sizeof(ui16);
            buf2 = swap_bytes_if_le(buf2);
            memcpy(dp, &buf2, 2); // not to violate strict aliasing
          }
        }
        else
        {
          ui32 buf4;
          char* dp = (char*)p->rec.marker_points;
          for (ui32 i = 0; i < p->rec.num_points; ++i, dp += 4)
          {
            result &= file->read(&buf4, sizeof(ui32)) == sizeof(ui32);
            buf4 = swap_bytes_if_le(buf4);
            memcpy(dp, &buf4, 4); // not to violate strict aliasing
          }
        }
        p->rec.prepare_for_decoding();
      }
 
      if (result == false)
        OJPH_ERROR(0x000501E3, "Error reading NLT marker, Location 3");
 
      return result;
    }

    const nlt_rec* param_nlt::get_nlt_rec(ui32 comp_num) const
    {
      assert(Cnlt == special_comp_num::ALL_COMPS);
      const param_nlt* p = get_nlt_object(comp_num);
      if (p && p->enabled)
        return &p->rec;
      else if (this->enabled)
        return &this->rec;
      else
        return NULL;
    }

    param_nlt* param_nlt::get_nlt_object(ui32 comp_num)
    {
      // cast object to constant
      const param_nlt* const_p = const_cast<const param_nlt*>(this);
      // call using the constant object, then cast to non-const
      return const_cast<param_nlt*>(const_p->get_nlt_object(comp_num));
    }

    const param_nlt* param_nlt::get_nlt_object(ui32 comp_num) const
    {
      // top level must be ALL_COMPS; cannot call this functions from
      // non-top nlt
      assert(Cnlt == special_comp_num::ALL_COMPS);
 
      const param_nlt* p = this;
      while (p && p->Cnlt != comp_num)
        p = p->next;
      return p;
    }

    param_nlt* param_nlt::add_object(ui32 comp_num)
    {
      assert(comp_num != special_comp_num::ALL_COMPS);
      assert(Cnlt == special_comp_num::ALL_COMPS);      // top level
      param_nlt* p = this;
      while (p->next != NULL) {
        assert(p->Cnlt != comp_num);
        p = p->next;
      }
      if (avail)
      {
        p->next = avail;
        avail = avail->next;
        p->next->init();
      }
      else
        p->next = new param_nlt;
      p = p->next;
      p->Cnlt = (ui16)comp_num;
      return p;
    }

    bool param_nlt::is_any_enabled() const
    {
      // check if any field is enabled
      const param_nlt* p = this;
      while (p && p->enabled == false)
        p = p->next;
      return (p != NULL);
    }

    void param_nlt::trim_non_existing_components(ui32 num_comps)
    {
      param_nlt* p = this->next;
      while (p) {
          if (p->enabled == true && p->Cnlt >= num_comps) {
            p->enabled = false;
            OJPH_INFO(0x00050161, "The NLT marker segment for the "
              "non-existing component %d has been removed.", p->Cnlt);
          }
        p = p->next;
      }
    }
 

    //
    //
    //
    //
    //

    bool param_sot::write(outfile_base *file, ui32 payload_len)
    {
      ui16 buf2;
      ui32 buf4;
      bool result = true;
 
      this->Psot = payload_len + 14; //inc. SOT marker, field & SOD
 
      buf2 = JP2K_MARKER::SOT;
      buf2 = swap_bytes_if_le(buf2);
      result &= file->write(&buf2, sizeof(ui16)) == sizeof(ui16);
      buf2 = swap_bytes_if_le(Lsot);
      result &= file->write(&buf2, sizeof(ui16)) == sizeof(ui16);
      buf2 = swap_bytes_if_le(Isot);
      result &= file->write(&buf2, sizeof(ui16)) == sizeof(ui16);
      buf4 = swap_bytes_if_le(Psot);
      result &= file->write(&buf4, sizeof(ui32)) == sizeof(ui32);
      result &= file->write(&TPsot, 1) == 1;
      result &= file->write(&TNsot, 1) == 1;
 
      return result;
    }

    bool param_sot::write(outfile_base *file, ui32 payload_len,
                          ui8 TPsot, ui8 TNsot)
    {
      ui32 buf4;
      ui16 buf2;
      bool result = true;
 
      buf2 = JP2K_MARKER::SOT;
      buf2 = swap_bytes_if_le(buf2);
      result &= file->write(&buf2, sizeof(ui16)) == sizeof(ui16);
      buf2 = swap_bytes_if_le(Lsot);
      result &= file->write(&buf2, sizeof(ui16)) == sizeof(ui16);
      buf2 = swap_bytes_if_le(Isot);
      result &= file->write(&buf2, sizeof(ui16)) == sizeof(ui16);
      buf4 = swap_bytes_if_le(payload_len + 14);
      result &= file->write(&buf4, sizeof(ui32)) == sizeof(ui32);
      result &= file->write(&TPsot, 1) == 1;
      result &= file->write(&TNsot, 1) == 1;
 
      return result;
    }

    bool param_sot::read(infile_base *file, bool resilient)
    {
      if (resilient)
      {
        if (file->read(&Lsot, 2) != 2)
        {
          OJPH_INFO(0x00050091, "error reading SOT marker");
          Lsot = 0; Isot = 0; Psot = 0; TPsot = 0; TNsot = 0;
          return false;
        }
        Lsot = swap_bytes_if_le(Lsot);
        if (Lsot != 10)
        {
          OJPH_INFO(0x00050092, "error in SOT length");
          Lsot = 0; Isot = 0; Psot = 0; TPsot = 0; TNsot = 0;
          return false;
        }
        if (file->read(&Isot, 2) != 2)
        {
          OJPH_INFO(0x00050093, "error reading tile index");
          Lsot = 0; Isot = 0; Psot = 0; TPsot = 0; TNsot = 0;
          return false;
        }
        Isot = swap_bytes_if_le(Isot);
        if (Isot == 0xFFFF)
        {
          OJPH_INFO(0x00050094, "tile index in SOT marker cannot be 0xFFFF");
          Lsot = 0; Isot = 0; Psot = 0; TPsot = 0; TNsot = 0;
          return false;
        }
        if (file->read(&Psot, 4) != 4)
        {
          OJPH_INFO(0x00050095, "error reading SOT marker");
          Lsot = 0; Isot = 0; Psot = 0; TPsot = 0; TNsot = 0;
          return false;
        }
        Psot = swap_bytes_if_le(Psot);
        if (file->read(&TPsot, 1) != 1)
        {
          OJPH_INFO(0x00050096, "error reading SOT marker");
          Lsot = 0; Isot = 0; Psot = 0; TPsot = 0; TNsot = 0;
          return false;
        }
        if (file->read(&TNsot, 1) != 1)
        {
          OJPH_INFO(0x00050097, "error reading SOT marker");
          Lsot = 0; Isot = 0; Psot = 0; TPsot = 0; TNsot = 0;
          return false;
        }
      }
      else
      {
        if (file->read(&Lsot, 2) != 2)
          OJPH_ERROR(0x00050091, "error reading SOT marker");
        Lsot = swap_bytes_if_le(Lsot);
        if (Lsot != 10)
          OJPH_ERROR(0x00050092, "error in SOT length");
        if (file->read(&Isot, 2) != 2)
          OJPH_ERROR(0x00050093, "error reading SOT tile index");
        Isot = swap_bytes_if_le(Isot);
        if (Isot == 0xFFFF)
          OJPH_ERROR(0x00050094, "tile index in SOT marker cannot be 0xFFFF");
        if (file->read(&Psot, 4) != 4)
          OJPH_ERROR(0x00050095, "error reading SOT marker");
        Psot = swap_bytes_if_le(Psot);
        if (file->read(&TPsot, 1) != 1)
          OJPH_ERROR(0x00050096, "error reading SOT marker");
        if (file->read(&TNsot, 1) != 1)
          OJPH_ERROR(0x00050097, "error reading SOT marker");
      }
      return true;
    }

    //
    //
    //
    //
    //

    void param_tlm::init(ui32 num_pairs, Ttlm_Ptlm_pair *store)
    {
      if (num_pairs > MAX_PAIRS_PER_SEG * MAX_SEGMENTS)
        OJPH_ERROR(0x000500B1, "Trying to store %d tileparts in TLM markers, "
                   "but at most %d can be indexed; a codestream can carry "
                   "%d TLM marker segments of %d entries each.",
                   num_pairs, MAX_PAIRS_PER_SEG * MAX_SEGMENTS,
                   MAX_SEGMENTS, MAX_PAIRS_PER_SEG);
      this->num_pairs = num_pairs;
      pairs = store;
    }

    void param_tlm::set_next_pair(ui16 Ttlm, ui32 Ptlm)
    {
      assert(next_pair_index < num_pairs);
      pairs[next_pair_index].Ttlm = Ttlm;
      pairs[next_pair_index].Ptlm = Ptlm + 14;
      ++next_pair_index;
    }

    bool param_tlm::write(outfile_base *file)
    {
      assert(next_pair_index == num_pairs);
      ui16 buf2;
      ui32 buf4;
      bool result = true;
 
      ui32 written = 0;
      ui32 z = 0;
      do
      {
        ui32 left = num_pairs - written;
        ui32 count = left < MAX_PAIRS_PER_SEG ? left : MAX_PAIRS_PER_SEG;
        ui8 Ztlm = (ui8)z;
        ui8 Stlm = 0x60; // 2-byte Ttlm, 4-byte Ptlm
 
        buf2 = JP2K_MARKER::TLM;
        buf2 = swap_bytes_if_le(buf2);
        result &= file->write(&buf2, sizeof(ui16)) == sizeof(ui16);
        buf2 = swap_bytes_if_le((ui16)(4 + 6 * count));
        result &= file->write(&buf2, sizeof(ui16)) == sizeof(ui16);
        result &= file->write(&Ztlm, 1) == 1;
        result &= file->write(&Stlm, 1) == 1;
        for (ui32 i = written; i < written + count; ++i)
        {
          buf2 = swap_bytes_if_le(pairs[i].Ttlm);
          result &= file->write(&buf2, sizeof(ui16)) == sizeof(ui16);
          buf4 = swap_bytes_if_le(pairs[i].Ptlm);
          result &= file->write(&buf4, sizeof(ui32)) == sizeof(ui32);
        }
        written += count;
        ++z;
      } while (written < num_pairs);
      return result;
    }

    //
    //
    //
    //
    //

    const param_dfs* param_dfs::get_dfs(int index) const
    {
      const param_dfs* p = this;
      while (p && p->Sdfs != index)
        p = p->next;
      return p;
    }

    param_dfs::dfs_dwt_type param_dfs::get_dwt_type(ui32 decomp_level) const
    {
      decomp_level = ojph_min(decomp_level, Ids);
      ui32 d = decomp_level - 1;          // decomp_level starts from 1
      ui32 idx = d >> 2;                  // complete bytes
      ui32 bits = d & 0x3;                // bit within the bytes
      ui32 val = (Ddfs[idx] >> (6 - 2 * bits)) & 0x3;
      return (dfs_dwt_type)val;
    }

    ui32 param_dfs::get_subband_idx(ui32 num_decompositions, ui32 resolution,
                                    ui32 subband) const
    {
      assert((resolution == 0 && subband == 0) ||
              (resolution > 0 && subband > 0 && subband < 4));
 
      ui32 ns[4] = { 0, 3, 1, 1 };
 
      ui32 idx = 0;
      if (resolution > 0)
      {
        idx = 0;
        ui32 i = 1;
        for (; i < resolution; ++i)
          idx += ns[get_dwt_type(num_decompositions - i + 1)];
        dfs_dwt_type t = get_dwt_type(num_decompositions - i + 1);
        idx += subband;
        if (t == VERT_DWT && subband == 2)
          --idx;
      }
 
      return idx;
    }

    point param_dfs::get_res_downsamp(ui32 skipped_resolutions) const
    {
      point factor(1, 1);
      ui32 decomp_level = 1;
      while (skipped_resolutions > 0)
      {
        param_dfs::dfs_dwt_type type = get_dwt_type(decomp_level);
        if (type == BIDIR_DWT)
        { factor.x *= 2; factor.y *= 2; }
        else if (type == HORZ_DWT)
          factor.x *= 2;
        else if (type == VERT_DWT)
          factor.y *= 2;
 
        ++decomp_level;
        --skipped_resolutions;
      }
      return factor;
    }

    bool param_dfs::read(infile_base *file)
    {
      if (Ldfs != 0) { // this param_dfs is used
        param_dfs* p = this;
        while (p->next != NULL)
          p = p->next;
        if (avail)
        {
          p->next = avail;
          avail = avail->next;
          p->next->init();
        }
        else
          p->next = new param_dfs;
        p = p->next;
        return p->read(file);
      }
 
      if (file->read(&Ldfs, 2) != 2)
        OJPH_ERROR(0x000500D1, "error reading DFS-Ldfs parameter");
      Ldfs = swap_bytes_if_le(Ldfs);
      if (file->read(&Sdfs, 2) != 2)
        OJPH_ERROR(0x000500D2, "error reading DFS-Sdfs parameter");
      Sdfs = swap_bytes_if_le(Sdfs);
      if (Sdfs > 15)
        OJPH_ERROR(0x000500D3, "The DFS-Sdfs parameter is %d, which is "
          "larger than the permissible 15", Sdfs);
      ui8 t, l_Ids = 0;
      if (file->read(&l_Ids, 1) != 1)
        OJPH_ERROR(0x000500D4, "error reading DFS-Ids parameter");
      if (l_Ids == 0)
        OJPH_ERROR(0x000500D8,
          "The value of the Ids member in the DFS marker segment cannot be 0");
      constexpr int max_Ddfs = sizeof(Ddfs) * 4;
      if (l_Ids > max_Ddfs)
        OJPH_INFO(0x000500D5, "The DFS-Ids parameter is %d; while this is "
          "valid, the number is unnessarily large -- you do not need more "
          "than %d.  Please contact me regarding this issue.",
          l_Ids, max_Ddfs);
      Ids = l_Ids < max_Ddfs ? l_Ids : max_Ddfs;
      for (int i = 0; i < Ids; i += 4)
        if (file->read(&Ddfs[i / 4], 1) != 1)
          OJPH_ERROR(0x000500D6, "error reading DFS-Ddfs parameters");
      for (int i = Ids; i < l_Ids; i += 4)
        if (file->read(&t, 1) != 1)
          OJPH_ERROR(0x000500D7, "error reading DFS-Ddfs parameters");
      return true;
    }

    //
    //
    //
    //
    //

    param_atk* param_atk::get_atk(int index)
    {
      assert(top_atk == NULL);
 
      if (Latk == 0)
      {
        // This atk object is not used, initialize it to either 0 (irv97)
        // or 1 (rev53), and use it.  If index is not 0 nor 1, then index
        // must have been read from file previously, otherwise it is an
        // error.
        if (index == 0) { this->init_irv97(); return this; }
        else if (index == 1) { this->init_rev53(); return this; }
      }
 
      param_atk* p = this;
      while (p && p->get_index() != index)
        p = p->next;
 
      if (p == NULL && (index == 0 || index == 1))
      {
        // The index was not found, add an atk object only if the index is
        // either 0 or 1
        p = add_object();
        if (index == 0)
          p->init_irv97();
        else if (index == 1)
          p->init_rev53();
      }
 
      return p;
    }

    bool param_atk::read_coefficient(infile_base *file, float &K, si32& bytes)
    {
      int coeff_type = get_coeff_type();
      if (coeff_type == 0) { // 8bit
        ui8 v;
        if (file->read(&v, 1) != 1) return false;
        bytes -= 1;
        K = v;
      }
      else if (coeff_type == 1) { // 16bit
        ui16 v;
        if (file->read(&v, 2) != 2) return false;
        bytes -= 2;
        K = swap_bytes_if_le(v);
      }
      else if (coeff_type == 2) { // float
        ui32 i;
        if (file->read(&i, sizeof(ui32)) != sizeof(ui32)) return false;
        bytes -= 4;
        i = swap_bytes_if_le(i);
        float f;
        memcpy(&f, &i, sizeof(float));
        K = f;
      }
      else if (coeff_type == 3) { // double
        ui64 i;
        if (file->read(&i, sizeof(ui64)) != sizeof(ui64)) return false;
        bytes -= 8;
        i = swap_bytes_if_le(i);
        double d;
        memcpy(&d, &i, sizeof(double));
        K = (float)d;
      }
      else if (coeff_type == 4) { // 128 bit float
        ui64 v, v1;
        if (file->read(&v, 8) != 8) return false;
        bytes -= 8;
        if (file->read(&v1, 8) != 8) return false; // v1 not needed
        bytes -= 8;
        v = swap_bytes_if_le(v);
 
        // convert the MSB of 128b float to 32b float
        // 32b float has 1 sign bit, 8 exponent (offset 127), 23 mantissa
        // 128b float has 1 sign bit, 15 exponent (offset 16383), 112 mantissa
        si32 e = (si32)((v >> 48) & 0x7FFF);   // exponent
        e -= 16383;
        e += 127;
        e = e & 0xFF;                          // removes MSBs if negative
        e <<= 23;                              // move bits to their location
        ui32 i = 0;
        i |= ((ui32)(v >> 32) & 0x80000000); // copy sign bit
        i |= (ui32)e;                        // copy exponent
        i |= (ui32)((v >> 25) & 0x007FFFFF); // copy 23 mantissa
        float f;
        memcpy(&f, &i, sizeof(float));
        K = f;
      }
      return true;
    }
 

    bool param_atk::read_coefficient(infile_base *file, si16 &K, si32& bytes)
    {
      int coeff_type = get_coeff_type();
      if (coeff_type == 0) {
        si8 v;
        if (file->read(&v, 1) != 1) return false;
        bytes -= 1;
        K = v;
      }
      else if (coeff_type == 1) {
        si16 v;
        if (file->read(&v, 2) != 2) return false;
        bytes -= 2;
        K = (si16)swap_bytes_if_le((ui16)v);
      }
      else
        return false;
      return true;
    }

    bool param_atk::read(infile_base *file)
    {
      if (Latk != 0) // this param_atk is used
        return add_object()->read(file);
 
      if (file->read(&Latk, 2) != 2)
        OJPH_ERROR(0x000500E1, "error reading ATK-Latk parameter");
      Latk = swap_bytes_if_le(Latk);
      si32 bytes = Latk - 2;
      ojph::ui16 temp_Satk;
      if (file->read(&temp_Satk, 2) != 2)
        OJPH_ERROR(0x000500E2, "error reading ATK-Satk parameter");
      bytes -= 2;
      temp_Satk = swap_bytes_if_le(temp_Satk);
      int tmp_idx = temp_Satk & 0xFF;
      if ((top_atk && top_atk->get_atk(tmp_idx) != NULL)
        || tmp_idx == 0 || tmp_idx == 1)
        OJPH_ERROR(0x000500F3, "ATK-Satk parameter sets ATK marker index to "
          "the illegal value of %d. ATK-Satk should be in (2-255) and, I "
          "believe, must not be repeated; otherwise, it would be unclear "
          "what marker segment must be employed when an index is repeated.",
          tmp_idx);
      Satk = temp_Satk;
      if (is_m_init0() == false)  // only even-indexed is supported
        OJPH_ERROR(0x000500E3, "ATK-Satk parameter sets m_init to 1, "
          "requiring odd-indexed subsequence in first reconstruction step, "
          "which is not supported yet.");
      if (is_whole_sample() == false)  // ARB filter not supported
        OJPH_ERROR(0x000500E4, "ATK-Satk parameter specified ARB filter, "
          "which is not supported yet.");
      if (is_reversible() && get_coeff_type() >= 2) // reversible & float
        OJPH_ERROR(0x000500E5, "ATK-Satk parameter does not make sense. "
          "It employs floats with reversible filtering.");
      if (is_using_ws_extension() == false)  // only sym. ext is supported
        OJPH_ERROR(0x000500E6, "ATK-Satk parameter requires constant "
          "boundary extension, which is not supported yet.");
      if (is_reversible() == false)
        if (read_coefficient(file, Katk, bytes) == false)
          OJPH_ERROR(0x000500E7, "error reading ATK-Katk parameter");
      if (file->read(&Natk, 1) != 1)
        OJPH_ERROR(0x000500E8, "error reading ATK-Natk parameter");
      bytes -= 1;
      if (Natk > max_steps) {
        if (d != d_store) // was this allocated -- very unlikely
          delete[] d;
        d = new lifting_step[Natk];
        max_steps = Natk;
      }
 
      if (is_reversible())
      {
        for (int s = 0; s < Natk; ++s)
        {
          if (file->read(&d[s].rev.Eatk, 1) != 1)
            OJPH_ERROR(0x000500E9, "error reading ATK-Eatk parameter");
          bytes -= 1;
          if (file->read(&d[s].rev.Batk, 2) != 2)
            OJPH_ERROR(0x000500EA, "error reading ATK-Batk parameter");
          bytes -= 2;
          d[s].rev.Batk = (si16)swap_bytes_if_le((ui16)d[s].rev.Batk);
          ui8 LCatk;
          if (file->read(&LCatk, 1) != 1)
            OJPH_ERROR(0x000500EB, "error reading ATK-LCatk parameter");
          bytes -= 1;
          if (LCatk == 0)
            OJPH_ERROR(0x000500EC, "Encountered a ATK-LCatk value of zero; "
              "something is wrong.");
          if (LCatk > 1)
            OJPH_ERROR(0x000500ED, "ATK-LCatk value greater than 1; "
              "that is, a multitap filter is not supported");
          if (read_coefficient(file, d[s].rev.Aatk, bytes) == false)
            OJPH_ERROR(0x000500EE, "Error reding ATK-Aatk parameter");
        }
      }
      else
      {
        for (int s = 0; s < Natk; ++s)
        {
          ui8 LCatk;
          if (file->read(&LCatk, 1) != 1)
            OJPH_ERROR(0x000500EF, "error reading ATK-LCatk parameter");
          bytes -= 1;
          if (LCatk == 0)
            OJPH_ERROR(0x000500F0, "Encountered a ATK-LCatk value of zero; "
              "something is wrong.");
          if (LCatk > 1)
            OJPH_ERROR(0x000500F1, "ATK-LCatk value greater than 1; "
              "that is, a multitap filter is not supported.");
          if (read_coefficient(file, d[s].irv.Aatk, bytes) == false)
            OJPH_ERROR(0x000500F2, "Error reding ATK-Aatk parameter");
        }
      }
      if (bytes != 0)
        OJPH_ERROR(0x000500F3, "The length of an ATK marker segment "
          "(ATK-Latk) is not correct");
 
      return true;
    }

    void param_atk::init_irv97()
    {
      Satk = 0x4a00;     // illegal because ATK = 0
      Katk = (float)1.230174104914001;
      Natk = 4;
      // next is (A-4) in T.801 second line
      Latk = (ui16)(5 + Natk + sizeof(float) * (1 + Natk));
      d[0].irv.Aatk = (float)0.443506852043971;
      d[1].irv.Aatk = (float)0.882911075530934;
      d[2].irv.Aatk = (float)-0.052980118572961;
      d[3].irv.Aatk = (float)-1.586134342059924;
    }

    void param_atk::init_rev53()
    {
      Satk = 0x5801;     // illegal because ATK = 1
      Natk = 2;
      // next is (A-4) in T.801 fourth line
      Latk = (ui16)(5 + 2 * Natk + sizeof(ui8) * (Natk + Natk));
      d[0].rev.Aatk = 1;
      d[0].rev.Batk = 2;
      d[0].rev.Eatk = 2;
      d[1].rev.Aatk = -1;
      d[1].rev.Batk = 1;
      d[1].rev.Eatk = 1;
    }

    param_atk* param_atk::add_object()
    {
      assert(top_atk == NULL);
      param_atk *p = this;
      while (p->next != NULL)
        p = p->next;
      if (avail)
      {
        p->next = avail;
        avail = avail->next;
      }
      else
        p->next = new param_atk;
      p = p->next;
      p->init(this);
      return p;
    }
 
  } // !local namespace
}  // !ojph namespace