diff --git a/src/apps/ojph_compress/ojph_compress.cpp b/src/apps/ojph_compress/ojph_compress.cpp index 05d77623..a5b026b6 100644 --- a/src/apps/ojph_compress/ojph_compress.cpp +++ b/src/apps/ojph_compress/ojph_compress.cpp @@ -701,8 +701,7 @@ int main(int argc, char * argv[]) { if (!reversible && quantization_step != -1.0f) codestream.access_qcd().set_irrev_quant(quantization_step); if (!reversible && qfactor != -1) - codestream.access_qcd().set_qfactor(0, - ojph::param_qcd::OJPH_COMP_Y, (ojph::ui8)qfactor); + codestream.access_qcd().set_qfactor((ojph::ui8)qfactor); if (profile_string[0] != '\0') codestream.set_profile(profile_string); codestream.set_tilepart_divisions(tileparts_at_resolutions, @@ -759,12 +758,8 @@ int main(int argc, char * argv[]) { cod.set_reversible(reversible); if (!reversible && quantization_step != -1.0f) codestream.access_qcd().set_irrev_quant(quantization_step); - if (!reversible && qfactor != -1) { - ojph::param_qcd qcd = codestream.access_qcd(); - qcd.set_qfactor(0, ojph::param_qcd::OJPH_COMP_Y, (ojph::ui8)qfactor); - qcd.set_qfactor(1, ojph::param_qcd::OJPH_COMP_CB, (ojph::ui8)qfactor); - qcd.set_qfactor(2, ojph::param_qcd::OJPH_COMP_CR, (ojph::ui8)qfactor); - } + if (!reversible && qfactor != -1) + codestream.access_qcd().set_qfactor((ojph::ui8)qfactor); codestream.set_planar(false); if (profile_string[0] != '\0') codestream.set_profile(profile_string); @@ -843,20 +838,9 @@ int main(int argc, char * argv[]) { cod.set_color_transform(employ_color_transform == 1); } cod.set_reversible(reversible); - if (!reversible && qfactor != -1) { - ojph::param_qcd qcd = codestream.access_qcd(); - if (num_comps == 1) { - qcd.set_qfactor(0, - ojph::param_qcd::OJPH_COMP_Y, (ojph::ui8)qfactor); - } else { - qcd.set_qfactor(0, - ojph::param_qcd::OJPH_COMP_Y, (ojph::ui8)qfactor); - qcd.set_qfactor(1, - ojph::param_qcd::OJPH_COMP_CB, (ojph::ui8)qfactor); - qcd.set_qfactor(2, - ojph::param_qcd::OJPH_COMP_CR, (ojph::ui8)qfactor); - } - } else if (!reversible) { + if (!reversible && qfactor != -1) + codestream.access_qcd().set_qfactor((ojph::ui8)qfactor); + else if (!reversible) { const float min_step = 1.0f / 16384.0f; if (quantization_step == -1.0f) quantization_step = min_step; @@ -939,18 +923,8 @@ int main(int argc, char * argv[]) { OJPH_ERROR(0x010000A3, "-qfactor is only supported for images with 1 or 3 " "components\n"); - ojph::param_qcd qcd = codestream.access_qcd(); - if (num_comps == 1) { - qcd.set_qfactor(0, - ojph::param_qcd::OJPH_COMP_Y, (ojph::ui8)qfactor); - } else { - qcd.set_qfactor(0, - ojph::param_qcd::OJPH_COMP_Y, (ojph::ui8)qfactor); - qcd.set_qfactor(1, - ojph::param_qcd::OJPH_COMP_CB, (ojph::ui8)qfactor); - qcd.set_qfactor(2, - ojph::param_qcd::OJPH_COMP_CR, (ojph::ui8)qfactor); - } + codestream.access_qcd().set_qfactor((ojph::ui8)qfactor); + } codestream.set_planar(false); if (profile_string[0] != '\0') @@ -1042,18 +1016,7 @@ int main(int argc, char * argv[]) { OJPH_ERROR(0x010000A4, "-qfactor is only supported for images with 1 or 3 " "components\n"); - ojph::param_qcd qcd = codestream.access_qcd(); - if (num_components == 1) { - qcd.set_qfactor(0, - ojph::param_qcd::OJPH_COMP_Y, (ojph::ui8)qfactor); - } else { - qcd.set_qfactor(0, - ojph::param_qcd::OJPH_COMP_Y, (ojph::ui8)qfactor); - qcd.set_qfactor(1, - ojph::param_qcd::OJPH_COMP_CB, (ojph::ui8)qfactor); - qcd.set_qfactor(2, - ojph::param_qcd::OJPH_COMP_CR, (ojph::ui8)qfactor); - } + codestream.access_qcd().set_qfactor((ojph::ui8)qfactor); } codestream.set_planar(true); if (profile_string[0] != '\0') @@ -1108,8 +1071,7 @@ int main(int argc, char * argv[]) { if (!reversible && quantization_step != -1.0f) codestream.access_qcd().set_irrev_quant(quantization_step); if (!reversible && qfactor != -1) - codestream.access_qcd().set_qfactor(0, - ojph::param_qcd::OJPH_COMP_Y, (ojph::ui8)qfactor); + codestream.access_qcd().set_qfactor((ojph::ui8)qfactor); codestream.set_planar(true); if (profile_string[0] != '\0') codestream.set_profile(profile_string); @@ -1155,18 +1117,7 @@ int main(int argc, char * argv[]) { OJPH_ERROR(0x010000A5, "-qfactor is only supported for images with 1 or 3 " "components\n"); - ojph::param_qcd qcd = codestream.access_qcd(); - if (num_comps == 1) { - qcd.set_qfactor(0, - ojph::param_qcd::OJPH_COMP_Y, (ojph::ui8)qfactor); - } else { - qcd.set_qfactor(0, - ojph::param_qcd::OJPH_COMP_Y, (ojph::ui8)qfactor); - qcd.set_qfactor(1, - ojph::param_qcd::OJPH_COMP_CB, (ojph::ui8)qfactor); - qcd.set_qfactor(2, - ojph::param_qcd::OJPH_COMP_CR, (ojph::ui8)qfactor); - } + codestream.access_qcd().set_qfactor((ojph::ui8)qfactor); } codestream.set_planar(false); if (profile_string[0] != '\0') diff --git a/src/core/codestream/ojph_params.cpp b/src/core/codestream/ojph_params.cpp index 2fa0024b..b87106df 100644 --- a/src/core/codestream/ojph_params.cpp +++ b/src/core/codestream/ojph_params.cpp @@ -43,10 +43,11 @@ #include "ojph_file.h" #include "ojph_params.h" -#include "ojph_visual_weighting.h" #include "ojph_params_local.h" #include "ojph_message.h" +#include "ojph_visual_weighting.h" + namespace ojph { //////////////////////////////////////////////////////////////////////////// @@ -414,6 +415,11 @@ namespace ojph { state->set_delta(delta); } + ////////////////////////////////////////////////////////////////////////// + void param_qcd::set_qfactor(ui8 qfactor) { + state->set_qfactor(qfactor); + } + ////////////////////////////////////////////////////////////////////////// void param_qcd::set_irrev_quant(ui32 comp_idx, float delta) { @@ -1156,12 +1162,13 @@ namespace ojph { ////////////////////////////////////////////////////////////////////////// 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) { @@ -1172,6 +1179,37 @@ namespace ojph { } } + // 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 @@ -1198,7 +1236,8 @@ namespace ojph { } ////////////////////////////////////////////////////////////////////////// - void param_qcd::make_quant_steps(ui32 comp_num, const param_cod &cod, const param_siz &siz) + void param_qcd::make_quant_steps(ui32 comp_num, const param_cod &cod, + const param_siz &siz) { if (this->is_init) OJPH_ERROR(0x00040001, "Quantization step sizes already initialized."); @@ -1223,12 +1262,16 @@ namespace ojph { ui32 t = ojph_min(16, bit_depth); this->base_delta = 1.0f / (float)(1 << t); } - this->set_irrev_quant(this->num_decomps); + if (qfactor == QFACTOR_UNSET) + this->set_irrev_quant(this->num_decomps); + else + this->set_qfactor_quant(this->num_decomps); } } ////////////////////////////////////////////////////////////////////////// - bool param_qcd::is_qcc_needed(ui32 comp_num, const param_cod &cod, const param_siz &siz) + bool param_qcd::is_qcc_needed(ui32 comp_num, const param_cod &cod, + const param_siz &siz) { if (! this->is_init) OJPH_ERROR(0x00040001, "Quantization step sizes not initialized."); @@ -1238,33 +1281,17 @@ namespace ojph { 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_delta(ui32 comp_idx, float delta) - { - assert(type == QCD_MAIN); - param_qcd *p = get_qcc(comp_idx); - if (p == NULL) - p = add_qcc_object(comp_idx); - p->set_delta(delta); } ////////////////////////////////////////////////////////////////////////// - void param_qcd::set_qfactor(ui32 comp_idx, ojph::param_qcd::comp_type ctype, ui8 qfactor) { - if (this->top_qcd != NULL) - OJPH_ERROR(0x00040401, "This method is not implemented for QCC."); + void param_qcd::set_qfactor(ui8 qfactor) { + assert(this->type == QCD_MAIN); if (qfactor < 1 || qfactor > 100) - OJPH_ERROR(0x00040403, "Qfactor must be between 1 and 100, but was set to %i.", qfactor); + OJPH_ERROR(0x00050181, "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; + this->qfactor = qfactor; } ////////////////////////////////////////////////////////////////////////// @@ -1320,39 +1347,66 @@ namespace ojph { int guard_bits = 1; Sqcd = (ui8)((guard_bits<<5)|0x2);//one guard bit, scalar quantization - // the following are used only when Qfactor is set - float qfactor_power = 0; - float delta_ref = 0; - float G_c = 1; - const open_htj2k::visual_weighting_params vp; - std::vector W_b; - if (this->qfactor != QFACTOR_UNSET) { - const open_htj2k::q_scaling qs = open_htj2k::q_to_delta(this->qfactor, (ui8) this->bit_depth); - qfactor_power = (float) qs.qfactor_power; - const open_htj2k::color_transform ct = open_htj2k::resolve_color_transform(vp, this->is_color_trans); - int comp_index = (int)this->ctype; - delta_ref = (float) (qs.delta_Q * open_htj2k::color_gain(ct, 0)); - G_c = (float) open_htj2k::color_gain(ct, comp_index); - - // chroma_format 0 = 4:4:4, 1 = 4:2:0, 2 = 4:2:2 - int sampling = 0; - if (this->sampling.x == 1 && this->sampling.y == 1) - { sampling = 0; } - else if (this->sampling.x == 2 && this->sampling.y == 2) - { sampling = 1; } - else if (this->sampling.x == 2 && this->sampling.y == 1) - { sampling = 2; } - else + // LL, HL, LH, HH, HL, LH, HH... + for (ui32 s = 0; s < (1 + num_decomps * 3); s++) + { + // compute square root of the enery gain factor W_g + float w_g = 1.0; + + if (num_decomps > 0) { - OJPH_ERROR(0x00050161, "Qfactor can only be used on components with 4:4:4, 4:2:2 or 4:2:0 sampling"); - } + //In C++, division result truncates towards zero + ui32 d = num_decomps - (ui32)(((int)s - 1) / 3); + float gain_l = sqrt_energy_gains::get_gain_l(d, false); + float gain_h = sqrt_energy_gains::get_gain_h(d - 1, false); - if (this->ctype == ojph::param_qcd::OJPH_COMP_Y) - W_b = open_htj2k::luma_visual_weights((ui8) num_decomps, vp); - else - W_b = open_htj2k::chroma_visual_weights((ui8) num_decomps, vp, comp_index, sampling, ct); + if (s == 0) + { w_g = gain_l * gain_l; } + else if ((s - 1) % 3 == 2) + { w_g = gain_h * gain_h; } + else + { w_g = gain_l * gain_h; } + } + float delta_b = base_delta / w_g; + encode_SPqcd(s, delta_b); } + } + + ////////////////////////////////////////////////////////////////////////// + void param_qcd::set_qfactor_quant(ui32 num_decomps) + { + int guard_bits = 1; + Sqcd = (ui8)((guard_bits<<5)|0x2);//one guard bit, scalar quantization + + // the following are used only when Qfactor is set + const open_htj2k::visual_weighting_params vp; + open_htj2k::color_transform ct = + open_htj2k::resolve_color_transform(vp, this->is_color_trans); + + const open_htj2k::q_scaling qs = open_htj2k::q_to_delta(this->qfactor, + (ui8)ojph_min(16, this->bit_depth)); + float qfactor_power = (float) qs.qfactor_power; + float delta_ref = (float) (qs.delta_Q * open_htj2k::color_gain(ct, 0)); + float G_c = (float) open_htj2k::color_gain(ct, this->ctype); + + int chroma_format = 0; + // chroma_format 0 = 4:4:4, 1 = 4:2:0, 2 = 4:2:2 + if (this->sampling.x == 1 && this->sampling.y == 1) + chroma_format = 0; + else if (this->sampling.x == 2 && this->sampling.y == 2) + chroma_format = 1; + else if (this->sampling.x == 2 && this->sampling.y == 1) + chroma_format = 2; + else + OJPH_ERROR(0x00050161, "Qfactor can only be used on components " + "with 4:4:4, 4:2:2 or 4:2:0 sampling"); + + open_htj2k::visual_weights weights = + open_htj2k::visual_weights::make_weights((ui8) num_decomps, vp, + comp_idx, chroma_format, ct); + size_t count = weights.get_count(); + const float* w = weights.get_weights(); // LL, HL, LH, HH, HL, LH, HH... for (ui32 s = 0; s < (1 + num_decomps * 3); s++) @@ -1376,27 +1430,27 @@ namespace ojph { } float delta_b; - if (this->qfactor == QFACTOR_UNSET) - { - delta_b = base_delta / w_g; - } - else - { - float w_b = (s == 0 || s > W_b.size()) ? (float) 1.0 : (float) pow(W_b[W_b.size() - s], qfactor_power); - delta_b = delta_ref / (w_g * w_b * G_c); - } - - int exp = 0, mantissa; - while (delta_b < 1.0f) - { exp++; delta_b *= 2.0f; } - mantissa = (int)round(delta_b * (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[s] = (ui16)((exp << 11) | mantissa); + float w_b = 1.0f; + if (s > 0 && s <= count) + w_b = std::pow(w[count - s], qfactor_power); + delta_b = delta_ref / (w_g * w_b * G_c); + encode_SPqcd(s, delta_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 { @@ -1793,6 +1847,33 @@ namespace ojph { 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 == NULL) + p = add_qcc_object(comp_idx); + p->set_delta(delta); + } + + ////////////////////////////////////////////////////////////////////////// + void param_qcd::set_qfactor(ui32 comp_idx, comp_type ctype, ui8 qfactor) + { + assert(this->type == QCD_MAIN); + + if (qfactor < 1 || qfactor > 100) + 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) { diff --git a/src/core/codestream/ojph_params_local.h b/src/core/codestream/ojph_params_local.h index 40f41173..ab25fa8b 100644 --- a/src/core/codestream/ojph_params_local.h +++ b/src/core/codestream/ojph_params_local.h @@ -706,6 +706,9 @@ namespace ojph { QFACTOR_UNSET = 0 }; + //////////////////////////////////////// + using comp_type = ojph::param_qcd::comp_type; + public: param_qcd(param_qcd* top_qcd = NULL, ui16 comp_idx = OJPH_QCD_DEFAULT) { avail = NULL; init(top_qcd, comp_idx); } @@ -722,11 +725,12 @@ namespace ojph { } void check_validity(const param_siz& siz, const param_cod& cod); - void make_quant_steps(ui32 comp_num, const param_cod &cod, const param_siz &siz); - bool is_qcc_needed(ui32 comp_num, const param_cod &cod, const param_siz &siz); + void make_quant_steps(ui32 comp_num, const param_cod &cod, + const param_siz &siz); + bool is_qcc_needed(ui32 comp_num, const param_cod &cod, + const param_siz &siz); void set_delta(float delta) { base_delta = delta; } - void set_delta(ui32 comp_idx, float delta); - void set_qfactor(ui32 comp_idx, ojph::param_qcd::comp_type ctype, ui8 qfactor); + void set_qfactor(ui8 qfactor); ui32 get_num_guard_bits() const; ui32 get_MAGB() const; ui32 get_Kmax(const param_dfs* dfs, ui32 num_decompositions, @@ -740,6 +744,8 @@ namespace ojph { void read(infile_base *file); void read_qcc(infile_base *file, ui32 num_comps); + void set_delta(ui32 comp_idx, float delta); + void set_qfactor(ui32 comp_idx, comp_type ctype, ui8 qfactor); param_qcd* get_qcc(ui32 comp_idx); const param_qcd* get_qcc(ui32 comp_idx) const; param_qcd* add_qcc_object(ui32 comp_idx); @@ -757,7 +763,7 @@ namespace ojph { num_subbands = 0; base_delta = -1.0f; qfactor = QFACTOR_UNSET; - ctype = ojph::param_qcd::OJPH_COMP_Y; + ctype = comp_type::OJPH_COMP_Y; sampling = ojph::point(1, 1); enabled = true; next = NULL; @@ -780,6 +786,8 @@ namespace ojph { void set_rev_quant(ui32 num_decomps, ui32 bit_depth, bool is_employing_color_transform); void set_irrev_quant(ui32 num_decomps); + void set_qfactor_quant(ui32 num_decomps); + void encode_SPqcd(ui32 subband_index, float delta); ui32 get_largest_Kmax() const; bool internal_write_qcc(outfile_base *file, ui32 num_comps); void trim_non_existing_components(ui32 num_comps); @@ -808,7 +816,7 @@ namespace ojph { float base_delta; // base quantization step size -- all other // step sizes are derived from it. ui8 qfactor; - ojph::param_qcd::comp_type ctype; + comp_type ctype; bool is_color_trans; ui32 num_decomps; ui32 bit_depth; diff --git a/src/core/codestream/ojph_visual_weighting.h b/src/core/codestream/ojph_visual_weighting.h index 8533b5b8..33c7870c 100644 --- a/src/core/codestream/ojph_visual_weighting.h +++ b/src/core/codestream/ojph_visual_weighting.h @@ -2,10 +2,10 @@ // All rights reserved. // // Redistribution and use in source and binary forms, with or without -// modification, are permitted provided that the following conditions are met: +// 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. +// 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 @@ -17,16 +17,18 @@ // // 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. +// 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. -#pragma once +#ifndef OJPH_VISUAL_WEIGHTS_H +#define OJPH_VISUAL_WEIGHTS_H // --------------------------------------------------------------------------- // Analytic per-subband visual (CSF) weighting for Qfactor quantization. @@ -42,29 +44,30 @@ // * mannos_sakrison reproduces W_b_Y to RMS ~0.029 at ref_ppd ~= 72 // * daly reproduces W_b_Y near the paper's ~1700 px viewing // distance but with a looser shape (RMS ~0.056) -// Header-only and dependency-free (only //) so it can -// be exercised by a standalone demo as well as by the marker code. +// Header-only and dependency-free (only //) +// so it can be exercised by a standalone demo as well as by the marker code. // --------------------------------------------------------------------------- +#include #include #include -#include namespace open_htj2k { enum class csf_model { - legacy_table, // exact Zeng Table 2 weights (default; bit-identical output) - mannos_sakrison, // analytic Mannos-Sakrison CSF - daly // analytic Daly CSF (light-adaptation form) + legacy_table, // exact Zeng Table 2 weights (default; bit-identical output) + mannos_sakrison, // analytic Mannos-Sakrison CSF + daly // analytic Daly CSF (light-adaptation form) }; // Color decorrelation actually in force, which drives BOTH the per-component -// synthesis gain and whether a QCC component is treated as chroma (opponent) or -// luma (e.g. undecorrelated RGB). Quantization error in component c is amplified -// by ||column_c(MCT^-1)||_2 in reconstructed (R,G,B) space, so the step is scaled -// by 1/G_c. With no MCT every component is independent (unit gain, luma role). +// synthesis gain and whether a QCC component is treated as chroma (opponent) +// or luma (e.g. undecorrelated RGB). Quantization error in component c is +// amplified by ||column_c(MCT^-1)||_2 in reconstructed (R,G,B) space, so the +// step is scaled by 1/G_c. With no MCT every component is independent +// (unit gain, luma role). enum class color_transform { - ict, // irreversible 9/7 MCT (YCbCr): ICT inverse column-norm gains, chroma roles + ict, // irrev. 9/7 MCT (YCbCr): ICT inverse column-norm gains, chroma roles none // no MCT: unit gains, luma role for every component (RGB / generic) }; @@ -74,16 +77,16 @@ struct visual_weighting_params { csf_model model = csf_model::legacy_table; // Reference pixels-per-degree at zoom 1.0. Calibrated so that // mannos_sakrison reproduces the legacy W_b_Y table (see header note). - double ref_ppd = 72.0; + float ref_ppd = 72.0f; // Display magnification. zoom > 1 (zoom-in) lowers the effective ppd, shifts // every subband to a lower cycles/degree, and flattens the weighting toward // 1.0 (i.e. toward flat MSE-optimal quantization) -- the correct limit when // the viewer is close / heavily magnified. - double zoom = 1.0; + float zoom = 1.0f; // Radial-frequency multiplier for the diagonal (HH) subband relative to the // horizontal/vertical (LH/HL) center. Geometric value is sqrt(2); the legacy // table behaves closer to ~1.25 (no oblique-effect penalty). - double hh_factor = 1.4142135623730951; + float hh_factor = std::sqrt(2.0f); }; // --- color synthesis gain -------------------------------------------------- @@ -95,18 +98,21 @@ struct visual_weighting_params { // historical 4-dp literals so the emitted QCD/QCC stays bit-identical. // Per-component synthesis gain. `none` (no MCT) -> unit gain (independent -// components). `ict` -> the inverse-ICT column norms {sqrt(3), 1.80511, 1.57340} -// rounded to 4 dp exactly as the legacy encoder stored them. -inline double color_gain(color_transform ct, int comp_index) { +// components). `ict` -> the inverse-ICT column norms +// {sqrt(3), 1.80511, 1.57340} rounded to 4 dp exactly as the legacy encoder +// stored them. +inline float color_gain(color_transform ct, int comp_index) { if (ct == color_transform::none) return 1.0; - static const double g[3] = {1.7321, 1.8051, 1.5734}; + static const float g[3] = {1.7321f, 1.8051f, 1.5734f}; return g[(comp_index >= 0 && comp_index < 3) ? comp_index : 0]; } // Resolve the effective transform for color scaling. Legacy mode reproduces the // historical encoder, which assumed ICT regardless of the real MCT; analytic // modes honor the transform actually applied (identity when the MCT is off). -inline color_transform resolve_color_transform(const visual_weighting_params &vp, bool mct_on) { +inline color_transform +resolve_color_transform(const visual_weighting_params &vp, bool mct_on) +{ if (vp.model == csf_model::legacy_table) return color_transform::ict; return mct_on ? color_transform::ict : color_transform::none; } @@ -118,73 +124,78 @@ inline color_transform resolve_color_transform(const visual_weighting_params &vp // estimate_qfactor so all three invert to identical steps (guarded by the // qfest_* round-trip tests). RI is the component dynamic range in bits. struct q_scaling { - double delta_Q; // reference step before basis / visual / color normalization - double qfactor_power; // exponent applied to each subband's visual weight + float delta_Q; // reference step before basis / visual / color normalization + float qfactor_power; // exponent applied to each subband's visual weight }; inline q_scaling q_to_delta(uint8_t qfactor, uint8_t RI) { const uint8_t t0 = 65, t1 = 97; - const double alpha_T0 = 0.04, alpha_T1 = 0.10; - const double M_T0 = 2.0 * (1.0 - t0 / 100.0); - const double M_T1 = 2.0 * (1.0 - t1 / 100.0); - const double M_Q = (qfactor < 50) ? 50.0 / qfactor : 2.0 * (1.0 - qfactor / 100.0); - double alpha_Q = alpha_T0; - double qfactor_power = 1.0; - if (qfactor >= t1) { - qfactor_power = 0.0; + const float alpha_T0 = 0.04f, alpha_T1 = 0.10f; + const float M_T0 = 2.0f * (1.0f - t0 / 100.0f); + const float M_T1 = 2.0f * (1.0f - t1 / 100.0f); + const float M_Q = + (qfactor < 50) ? 50.0f / qfactor : 2.0f * (1.0f - qfactor / 100.0f); + float alpha_Q = alpha_T0; + float qfactor_power = 1.0; + if (qfactor >= t1) + { + qfactor_power = 0.0f; alpha_Q = alpha_T1; - } else if (qfactor > t0) { - qfactor_power = (std::log(M_T1) - std::log(M_Q)) / (std::log(M_T1) - std::log(M_T0)); + } + else if (qfactor > t0) + { + qfactor_power = + (std::log(M_T1) - std::log(M_Q)) / (std::log(M_T1) - std::log(M_T0)); alpha_Q = alpha_T1 * std::pow(alpha_T0 / alpha_T1, qfactor_power); } // eps0 = sqrt(1/2) / 2^RI. Use ldexp rather than `1 << RI`: bit-identical for // the encoder's RI (<= 16 bpp) but well-defined for any RI, so a crafted // high-bit-depth file fed to estimate_qfactor cannot trigger signed-shift UB. - const double eps0 = std::sqrt(0.5) * std::ldexp(1.0, -static_cast(RI)); + const float eps0 = std::sqrt(0.5f) * std::ldexp(1.0f, -static_cast(RI)); return {alpha_Q * M_Q + eps0, qfactor_power}; } // --- CSF models (spatial frequency f in cycles/degree) --------------------- -inline double mannos_sakrison_csf(double f) { - return 2.6 * (0.0192 + 0.114 * f) * std::exp(-std::pow(0.114 * f, 1.1)); +inline float mannos_sakrison_csf(float f) { + return 2.6f * (0.0192f + 0.114f * f) * std::exp(-std::pow(0.114f * f, 1.1f)); } -inline double daly_csf(double f) { +inline float daly_csf(float f) { // Daly 1993 CSF core, light-adaptation form at L = 100 cd/m^2. The low- // frequency image-area term is omitted because the low side is clamped flat. - constexpr double L = 100.0, eps = 0.9; - const double A = 0.801 * std::pow(1.0 + 0.7 / L, -0.2); - const double B = 0.3 * std::pow(1.0 + 100.0 / L, 0.15); - const double x = B * eps * f; - return A * eps * f * std::exp(-x) * std::sqrt(1.0 + 0.06 * std::exp(x)); + constexpr float L = 100.0f, eps = 0.9f; + const float A = 0.801f * std::pow(1.0f + 0.7f / L, -0.2f); + const float B = 0.3f * std::pow(1.0f + 100.0f / L, 0.15f); + const float x = B * eps * f; + return A * eps * f * std::exp(-x) * std::sqrt(1.0f + 0.06f * std::exp(x)); } -inline double csf_value(double f, csf_model m) { +inline float csf_value(float f, csf_model m) { return (m == csf_model::daly) ? daly_csf(f) : mannos_sakrison_csf(f); } -// Guard a visual weight against extreme viewing conditions (very large effective -// ppd) where the CSF underflows to 0 or Daly's exp() overflows to NaN. Either -// would make the step `delta/(basis*w*G)` divide by zero / propagate NaN and -// emit a corrupt QCD/QCC. Replaces only non-finite or non-positive values, so it -// is a no-op for every realistic weight (bit-identical normal output). -inline double finite_weight(double w) { - return (std::isfinite(w) && w > 0.0) ? w : 1e-300; -} +// Guard a visual weight against extreme viewing conditions (very large +// effective ppd) where the CSF underflows to 0 or Daly's exp() overflows to +// NaN. Either would make the step `delta/(basis*w*G)` divide by +// zero / propagate NaN and emit a corrupt QCD/QCC. Replaces only non-finite +// or non-positive values, so it is a no-op for every realistic weight +// (bit-identical normal output). +inline float finite_weight(float w) +{ return (std::isfinite(w) && w > 0.0f) ? w : 1e-38f; } // Peak (frequency, amplitude) of a CSF, found once by a coarse scan. Used to // normalize the weight to <= 1 and to clamp the band-pass low side to flat. struct csf_peak_t { - double f_peak; - double h_peak; + float f_peak; + float h_peak; }; inline csf_peak_t csf_peak(csf_model m) { auto scan = [](csf_model mm) { - csf_peak_t pk{0.0, 0.0}; + csf_peak_t pk{0.0f, 0.0f}; for (int i = 1; i < 100000; ++i) { - const double f = i * 0.002; // 0.002 .. 200 cpd - const double v = csf_value(f, mm); + const float f = (float)i * 0.002f; // 0.002 .. 200 cpd + const float v = csf_value(f, mm); if (v > pk.h_peak) { pk.h_peak = v; pk.f_peak = f; @@ -192,145 +203,213 @@ inline csf_peak_t csf_peak(csf_model m) { } return pk; }; - // Peak depends only on the model -> compute each once (thread-safe magic statics). + // Peak depends only on the model -> compute each once + // (thread-safe magic statics). static const csf_peak_t mannos_pk = scan(csf_model::mannos_sakrison); static const csf_peak_t daly_pk = scan(csf_model::daly); return (m == csf_model::daly) ? daly_pk : mannos_pk; } // Normalized CSF weight: flat (= 1) at/below the peak, falling CSF above it. -inline double csf_weight(double f, csf_model m, const csf_peak_t &pk) { - if (f <= pk.f_peak) return 1.0; +inline float csf_weight(float f, csf_model m, const csf_peak_t &pk) { + if (f <= pk.f_peak) return 1.0f; return finite_weight(csf_value(f, m) / pk.h_peak); } -// Square-root-domain luminance visual weights, one per detail subband, in the -// QCD/wmse build order (per level, finest first): [HH_l, LH_l, HL_l]. Length -// is 3 * dwt_levels for analytic models; the LL band is handled by the caller -// (weight 1.0). In legacy mode the historical 15-entry table is returned -// verbatim, so the caller's existing out-of-range guard still applies. -inline std::vector luma_visual_weights(uint8_t dwt_levels, - const visual_weighting_params &vp) { - if (vp.model == csf_model::legacy_table) { - // Zeng et al., Table 2, Y column (= sqrt of the MSE-domain weights). - return {0.0901, 0.2758, 0.2758, 0.7018, 0.8378, 0.8378, 1.0000, 1.0000, - 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000}; - } - - const csf_peak_t pk = csf_peak(vp.model); - const double zoom = (vp.zoom > 0.0) ? vp.zoom : 1.0; - const double ppd = vp.ref_ppd / zoom; // zoom-in lowers effective ppd - const double f_N = ppd / 2.0; // Nyquist in cycles/degree - - std::vector w; - w.reserve(static_cast(3) * dwt_levels); - for (uint8_t lvl = 1; lvl <= dwt_levels; ++lvl) { - // Geometric-mean radial center of octave band [f_N/2^lvl, f_N/2^(lvl-1)]. - const double f_r = f_N * std::pow(2.0, -static_cast(lvl)) * std::sqrt(2.0); - const double f_hh = f_r * vp.hh_factor; - w.push_back(csf_weight(f_hh, vp.model, pk)); // HH - w.push_back(csf_weight(f_r, vp.model, pk)); // LH - w.push_back(csf_weight(f_r, vp.model, pk)); // HL - } - return w; -} - // --- chroma (chrominance) analytic weighting ------------------------------- // -// Chrominance CSF is low-pass (not band-pass like luminance) and rolls off well -// below the luminance CSF. A single per-channel CSF, evaluated at each subband's -// per-axis radial frequency, reproduces ALL of the historical 4:4:4 / 4:2:0 / -// 4:2:2 QCC tables: chroma subsampling is just a frequency shift (sx, sy) -// applied to the same CSF, which also yields the 4:2:2 LH != HL anisotropy for -// free. Parameters below were fit to the 4:4:4 Cb/Cr tables at ref_ppd = 72 -// (see scripts/csf_fit_chroma.py): Cb RMS 0.010, Cr RMS 0.021; predicted 4:2:0 / 4:2:2 -// RMS <= 0.064. The chroma CSF is independent of vp.model's luminance choice. +// Chrominance CSF is low-pass (not band-pass like luminance) and rolls off +// well below the luminance CSF. A single per-channel CSF, evaluated at each +// subband's per-axis radial frequency, reproduces ALL of the historical +// 4:4:4 / 4:2:0 / 4:2:2 QCC tables: chroma subsampling is just a frequency +// shift (sx, sy) applied to the same CSF, which also yields the +// 4:2:2 LH != HL anisotropy for free. Parameters below were fit to the +// 4:4:4 Cb/Cr tables at ref_ppd = 72 (see scripts/csf_fit_chroma.py): +// Cb RMS 0.010, Cr RMS 0.021; predicted 4:2:0 / 4:2:2 RMS <= 0.064. +// The chroma CSF is independent of vp.model's luminance choice. // Stretched-exponential low-pass chroma CSF (= 1 at DC). f in cycles/degree. -inline double chroma_csf(double f, double a, double b) { return std::exp(-std::pow(a * f, b)); } +inline float chroma_csf(float f, float a, float b) +{ return std::exp(-std::pow(a * f, b)); } // Calibrated chroma CSF parameters per opponent channel (at ref_ppd = 72). struct chroma_csf_params { - double a, b; + float a, b; }; inline chroma_csf_params chroma_params_for(int comp_index) { - // comp_index 1 = Cb (blue-yellow, steeper roll-off); otherwise Cr (red-green). - return (comp_index == 1) ? chroma_csf_params{0.1173, 0.840} : chroma_csf_params{0.0699, 1.050}; + // comp_index 1 = Cb (blue-yellow, steeper roll-off); + // otherwise Cr (red-green). + return (comp_index == 1) ? + chroma_csf_params{0.1173f, 0.840f} : chroma_csf_params{0.0699f, 1.050f}; } -// Exact historical QCC table row (sqrt-domain). comp_index 1 = Cb, else Cr; -// chroma_format 0 = 4:4:4, 1 = 4:2:0, 2 = 4:2:2 (matches j2kmarkers YCC*). -inline std::vector legacy_chroma_row(int comp_index, int chroma_format) { - static const double Cb444[15] = {0.0263, 0.0863, 0.0863, 0.1362, 0.2564, 0.2564, 0.3346, 0.4691, - 0.4691, 0.5444, 0.6523, 0.6523, 0.7078, 0.7797, 0.7797}; - static const double Cr444[15] = {0.0773, 0.1835, 0.1835, 0.2598, 0.4130, 0.4130, 0.5040, 0.6464, - 0.6464, 0.7220, 0.8254, 0.8254, 0.8769, 0.9424, 0.9424}; - static const double Cb420[15] = {0.1362, 0.2564, 0.2564, 0.3346, 0.4691, 0.4691, 0.5444, 0.6523, - 0.6523, 0.7078, 0.7797, 0.7797, 1.0000, 1.0000, 1.0000}; - static const double Cr420[15] = {0.2598, 0.4130, 0.4130, 0.5040, 0.6464, 0.6464, 0.7220, 0.8254, - 0.8254, 0.8769, 0.9424, 0.9424, 1.0000, 1.0000, 1.0000}; - static const double Cb422[15] = {0.0863, 0.0863, 0.2564, 0.2564, 0.2564, 0.4691, 0.4691, 0.4691, - 0.6523, 0.6523, 0.6523, 0.7797, 0.7797, 0.7797, 1.0000}; - static const double Cr422[15] = {0.1835, 0.1835, 0.4130, 0.4130, 0.4130, 0.6464, 0.6464, 0.6464, - 0.8254, 0.8254, 0.8254, 0.9424, 0.9424, 0.9424, 1.0000}; - const double *r; - switch (chroma_format) { - case 1: r = (comp_index == 1) ? Cb420 : Cr420; break; // 4:2:0 - case 2: r = (comp_index == 1) ? Cb422 : Cr422; break; // 4:2:2 - default: r = (comp_index == 1) ? Cb444 : Cr444; break; // 4:4:4 - } - return std::vector(r, r + 15); -} +// structure for our visual weights +class visual_weights { -// Square-root-domain chroma visual weights, one per detail subband, in QCC order -// (per level, finest first): [HH_l, LH_l, HL_l]. legacy_table returns the -// historical row verbatim (bit-identical). Analytic models fold chroma -// subsampling into the effective horizontal/vertical ppd, so LH (vertical -// detail) and HL (horizontal detail) diverge under 4:2:2 as they should. -inline std::vector chroma_visual_weights(uint8_t dwt_levels, - const visual_weighting_params &vp, int comp_index, - int chroma_format, - color_transform ct = color_transform::ict) { - if (vp.model == csf_model::legacy_table) { - return legacy_chroma_row(comp_index, chroma_format); - } - // Without a luma/chroma decorrelating transform this component carries - // luminance (e.g. a raw RGB channel), so it takes the luminance CSF, never - // the chroma roll-off. Such components are not subsampled (full resolution). - if (ct == color_transform::none) { - return luma_visual_weights(dwt_levels, vp); - } + public: + static visual_weights make_weights(uint8_t dwt_levels, + const visual_weighting_params &vp, int comp_index, int chroma_format = 0, + color_transform ct = color_transform::ict) + { + if (comp_index == 0) + return make_luma_weights(dwt_levels, vp); + else + return + make_chroma_weights(dwt_levels, vp, comp_index, chroma_format, ct); + } - double sx = 1.0, sy = 1.0; // horizontal/vertical chroma subsampling factors - if (chroma_format == 1) { // 4:2:0 - sx = 2.0; - sy = 2.0; - } else if (chroma_format == 2) { // 4:2:2 - sx = 2.0; - sy = 1.0; - } + const float* get_weights() const { return weights; } + size_t get_count() const { return count; } - const chroma_csf_params cp = chroma_params_for(comp_index); - const double zoom = (vp.zoom > 0.0) ? vp.zoom : 1.0; - const double f_N = (vp.ref_ppd / zoom) / 2.0; // luma Nyquist (cpd) - const double f_Nx = f_N / sx; // chroma horizontal Nyquist - const double f_Ny = f_N / sy; // chroma vertical Nyquist - - std::vector w; - w.reserve(static_cast(3) * dwt_levels); - for (uint8_t lvl = 1; lvl <= dwt_levels; ++lvl) { - const double dx = f_Nx * std::pow(2.0, -static_cast(lvl)) * std::sqrt(2.0); - const double dy = f_Ny * std::pow(2.0, -static_cast(lvl)) * std::sqrt(2.0); - // HH radial = geometric per-axis diagonal, scaled by hh_factor relative to the - // isotropic sqrt(2) so the same knob tunes luma and chroma (default sqrt(2) = no change). - const double f_hh = (vp.hh_factor / std::sqrt(2.0)) * std::sqrt(dx * dx + dy * dy); - // chroma_csf = exp(-(a*f)^b) is <= 1 by construction; finite_weight only - // guards the underflow-to-0 case at extreme effective ppd. - w.push_back(finite_weight(chroma_csf(f_hh, cp.a, cp.b))); // HH - w.push_back(finite_weight(chroma_csf(dy, cp.a, cp.b))); // LH (vertical detail) - w.push_back(finite_weight(chroma_csf(dx, cp.a, cp.b))); // HL (horizontal detail) - } - return w; -} + private: // Either legacy or analytic weights + // Square-root-domain luminance visual weights, one per detail subband, in + // the QCD/wmse build order (per level, finest first): [HH_l, LH_l, HL_l]. + // Length is 3 * dwt_levels for analytic models; the LL band is handled by + // the caller (weight 1.0). In legacy mode the historical 15-entry table is + // returned verbatim, so the caller's existing out-of-range guard still + // applies. + static visual_weights make_luma_weights(uint8_t dwt_levels, + const visual_weighting_params &vp) + { + visual_weights w; + if (vp.model == csf_model::legacy_table) { + // Zeng et al., Table 2, Y column (= sqrt of the MSE-domain weights). + float a[15] = {0.0901f, 0.2758f, 0.2758f, 0.7018f, 0.8378f, 0.8378f, + 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f}; + w.set_weights(a, sizeof(a)/sizeof(float)); + return w; + } + + const csf_peak_t pk = csf_peak(vp.model); + const float zoom = (vp.zoom > 0.0f) ? vp.zoom : 1.0f; + const float ppd = vp.ref_ppd / zoom; // zoom-in lowers effective ppd + const float f_N = ppd / 2.0f; // Nyquist in cycles/degree + + for (uint8_t lvl = 1; lvl <= dwt_levels; ++lvl) { + //Geometric-mean radial center of octave band [f_N/2^lvl,f_N/2^(lvl-1)] + const float f_r = f_N * std::pow(2.0f, -(float)(lvl)) * (float)M_SQRT2; + const float f_hh = f_r * vp.hh_factor; + w.push_back(csf_weight(f_hh, vp.model, pk)); // HH + w.push_back(csf_weight(f_r, vp.model, pk)); // LH + w.push_back(csf_weight(f_r, vp.model, pk)); // HL + } + return w; + } + + // Square-root-domain chroma visual weights, one per detail subband, in QCC + // order (per level, finest first): [HH_l, LH_l, HL_l]. legacy_table returns + // the historical row verbatim (bit-identical). Analytic models fold chroma + // subsampling into the effective horizontal/vertical ppd, so LH (vertical + // detail) and HL (horizontal detail) diverge under 4:2:2 as they should. + static visual_weights make_chroma_weights(uint8_t dwt_levels, + const visual_weighting_params &vp, int comp_index, int chroma_format, + color_transform ct = color_transform::ict) + { + if (vp.model == csf_model::legacy_table) { + return legacy_chroma_row(comp_index, chroma_format); + } + // Without a luma/chroma decorrelating transform this component carries + // luminance (e.g. a raw RGB channel), so it takes the luminance CSF, + // never the chroma roll-off. Such components are not subsampled + // (full resolution). + if (ct == color_transform::none) { + return make_luma_weights(dwt_levels, vp); + } + + float sx = 1.0f, sy = 1.0f; // horz/vert chroma subsampling factors + if (chroma_format == 1) { // 4:2:0 + sx = 2.0f; + sy = 2.0f; + } else if (chroma_format == 2) { // 4:2:2 + sx = 2.0f; + sy = 1.0f; + } + + const chroma_csf_params cp = chroma_params_for(comp_index); + const float zoom = (vp.zoom > 0.0f) ? vp.zoom : 1.0f; + const float f_N = (vp.ref_ppd / zoom) / 2.0f; // luma Nyquist (cpd) + const float f_Nx = f_N / sx; // chroma horizontal Nyquist + const float f_Ny = f_N / sy; // chroma vertical Nyquist + + visual_weights w; + for (uint8_t lvl = 1; lvl <= dwt_levels; ++lvl) { + const float dx = f_Nx * std::pow(2.0f, -(float)(lvl)) * (float)M_SQRT2; + const float dy = f_Ny * std::pow(2.0f, -(float)(lvl)) * (float)M_SQRT2; + // HH radial = geometric per-axis diagonal, scaled by hh_factor relative + // to the isotropic sqrt(2) so the same knob tunes luma and chroma + // (default sqrt(2) = no change). + const float f_hh = + (vp.hh_factor / std::sqrt(2.0f)) * std::sqrt(dx * dx + dy * dy); + // chroma_csf = exp(-(a*f)^b) is <= 1 by construction; finite_weight + // only guards the underflow-to-0 case at extreme effective ppd. + w.push_back(finite_weight(chroma_csf(f_hh, cp.a, cp.b))); // HH + w.push_back(finite_weight(chroma_csf(dy, cp.a, cp.b))); // LH + w.push_back(finite_weight(chroma_csf(dx, cp.a, cp.b))); // HL + } + return w; + } + + private: // Legacy chroma + // Exact historical QCC table row (sqrt-domain). comp_index 1 = Cb, else Cr; + // chroma_format 0 = 4:4:4, 1 = 4:2:0, 2 = 4:2:2 (matches j2kmarkers YCC*). + static visual_weights legacy_chroma_row(int comp_index, int chroma_format) + { + static const float Cb444[15] = + {0.0263f, 0.0863f, 0.0863f, 0.1362f, 0.2564f, 0.2564f, 0.3346f, 0.4691f, + 0.4691f, 0.5444f, 0.6523f, 0.6523f, 0.7078f, 0.7797f, 0.7797f}; + static const float Cr444[15] = + {0.0773f, 0.1835f, 0.1835f, 0.2598f, 0.4130f, 0.4130f, 0.5040f, 0.6464f, + 0.6464f, 0.7220f, 0.8254f, 0.8254f, 0.8769f, 0.9424f, 0.9424f}; + static const float Cb420[15] = + {0.1362f, 0.2564f, 0.2564f, 0.3346f, 0.4691f, 0.4691f, 0.5444f, 0.6523f, + 0.6523f, 0.7078f, 0.7797f, 0.7797f, 1.0f, 1.0f, 1.0f}; + static const float Cr420[15] = + {0.2598f, 0.4130f, 0.4130f, 0.5040f, 0.6464f, 0.6464f, 0.7220f, 0.8254f, + 0.8254f, 0.8769f, 0.9424f, 0.9424f, 1.0f, 1.0f, 1.0f}; + static const float Cb422[15] = + {0.0863f, 0.0863f, 0.2564f, 0.2564f, 0.2564f, 0.4691f, 0.4691f, 0.4691f, + 0.6523f, 0.6523f, 0.6523f, 0.7797f, 0.7797f, 0.7797f, 1.0f}; + static const float Cr422[15] = + {0.1835f, 0.1835f, 0.4130f, 0.4130f, 0.4130f, 0.6464f, 0.6464f, 0.6464f, + 0.8254f, 0.8254f, 0.8254f, 0.9424f, 0.9424f, 0.9424f, 1.0f}; + const float *r; + switch (chroma_format) { + case 1: r = (comp_index == 1) ? Cb420 : Cr420; break; // 4:2:0 + case 2: r = (comp_index == 1) ? Cb422 : Cr422; break; // 4:2:2 + default: r = (comp_index == 1) ? Cb444 : Cr444; break; // 4:4:4 + } + visual_weights result; + result.set_weights(r, 15); + return result; + } + + private: + visual_weights() : weights{}, count(0) {} + + void push_back(float val) { + assert(count < max_weight_size); + weights[count++] = val; + } + + void set_weights(const float* src, std::size_t count) + { + assert(count <= max_weight_size); + this->count = count; + for (std::size_t i = 0; i < count; ++i) + weights[i] = src[i]; + } + + private: + // Maximum number of DWT levels supported. + constexpr static uint8_t max_dwt_levels = 32; + // Maximum number of visual weight entries (3 subbands per level). + constexpr static size_t max_weight_size = 3 * max_dwt_levels; + + private: + float weights[max_weight_size]; + std::size_t count; +}; } // namespace open_htj2k + +#endif // !OJPH_VISUAL_WEIGHTS_H diff --git a/src/core/openjph/ojph_params.h b/src/core/openjph/ojph_params.h index 4c7dafbf..4d0f284b 100644 --- a/src/core/openjph/ojph_params.h +++ b/src/core/openjph/ojph_params.h @@ -169,11 +169,19 @@ namespace ojph { class OJPH_EXPORT param_qcd { public: - enum comp_type : ui8 { - OJPH_COMP_Y, - OJPH_COMP_CB, - OJPH_COMP_CR + enum comp_type : ui8 { // Note the numbers are used by the code + OJPH_COMP_Y = 0, + OJPH_COMP_CB = 1, + OJPH_COMP_CR = 2, + OJPH_COMP_UNDEFINED = 0xFF }; + static comp_type ui8_2_comp_type(ui8 c) + { + if (c >= OJPH_COMP_Y && c <= OJPH_COMP_CR) + return static_cast(c); + else + return OJPH_COMP_UNDEFINED; + } param_qcd(local::param_qcd* p) : state(p) {} @@ -187,6 +195,17 @@ namespace ojph { */ void set_irrev_quant(float delta); + /** + * @brief Sets Qfactor + * + * Setting Qfactor takes precedence over setting an irreversible + * quantization base delta + * + * @param qfactor Compression quality as an integer between + * 1 (worst quality) and 100 (best quality) + */ + void set_qfactor(ui8 qfactor); + /** * @brief Set the irreversible quantization base delta for a specific * component @@ -203,11 +222,14 @@ namespace ojph { /** * @brief Sets Qfactor for a specific component. * - * Setting Qfactor takes precedence over setting an irreversible quantization base delta + * Setting Qfactor takes precedence over setting an irreversible + * quantization base delta * * @param comp_idx Component index - * @param ctype Indicates whether the component is a Y, Cb or Cr channel, after the ICT if present - * @param qfactor Compression quality as an integer between 1 (worst quality) and 100 (best quality) + * @param ctype Indicates whether the component is a Y, Cb or Cr channel, + * after the ICT if present + * @param qfactor Compression quality as an integer between + * 1 (worst quality) and 100 (best quality) */ void set_qfactor(ui32 comp_idx, comp_type ctype, ui8 qfactor);