From 568808aadb789b870879daf8e8668792738569e0 Mon Sep 17 00:00:00 2001 From: Aous Naman Date: Wed, 1 Jul 2026 20:05:31 +1000 Subject: [PATCH 1/8] For visual weighting, this replaces vector with array, to give rid of the new/delete. A small change for bit-depth-based quantization; now it is limited to 2^-16, if not explicitly chosen. --- src/core/codestream/ojph_params.cpp | 27 +- src/core/codestream/ojph_visual_weighting.h | 324 ++++++++++++-------- 2 files changed, 206 insertions(+), 145 deletions(-) diff --git a/src/core/codestream/ojph_params.cpp b/src/core/codestream/ojph_params.cpp index 886b4b6d..cc9d8c75 100644 --- a/src/core/codestream/ojph_params.cpp +++ b/src/core/codestream/ojph_params.cpp @@ -1319,11 +1319,13 @@ namespace ojph { float delta_ref = 0; float G_c = 1; const open_htj2k::visual_weighting_params vp; - std::vector W_b; + open_htj2k::visual_weights W_b; if (this->qfactor != QFACTOR_UNSET) { - const open_htj2k::q_scaling qs = open_htj2k::q_to_delta(this->qfactor, (ui8) this->bit_depth); + const open_htj2k::q_scaling qs = open_htj2k::q_to_delta(this->qfactor, + (ui8)ojph_min(16, 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); + 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); @@ -1331,20 +1333,20 @@ namespace ojph { // 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; } + sampling = 0; else if (this->sampling.x == 2 && this->sampling.y == 2) - { sampling = 1; } + sampling = 1; else if (this->sampling.x == 2 && this->sampling.y == 1) - { sampling = 2; } + sampling = 2; else - { - OJPH_ERROR(0x00050161, "Qfactor can only be used on components with 4:4:4, 4:2:2 or 4:2:0 sampling"); - } + 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 == 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); + W_b = open_htj2k::chroma_visual_weights((ui8) num_decomps, vp, + comp_index, sampling, ct); } @@ -1371,12 +1373,11 @@ 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); + float w_b = (s == 0 || s > W_b.size()) ? + 1.0f : pow(W_b[W_b.size() - s], qfactor_power); delta_b = delta_ref / (w_g * w_b * G_c); } diff --git a/src/core/codestream/ojph_visual_weighting.h b/src/core/codestream/ojph_visual_weighting.h index 8533b5b8..2d288a7a 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,14 +17,15 @@ // // 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 @@ -42,29 +43,65 @@ // * 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 +#include namespace open_htj2k { +// Maximum number of DWT levels supported. +constexpr uint8_t max_dwt_levels = 32; +// Maximum number of visual weight entries (3 subbands per level). +constexpr size_t max_weight_size = 3 * max_dwt_levels; +// visual array type +using visual_weighting_array = std::array; +// structure for our visual weights +class visual_weights { + public: + visual_weights() : weights(), count(0) {} + + std::size_t size() const { return count; } + float operator[](std::size_t idx) { return weights[idx]; } + const float operator[](std::size_t idx) const { return weights[idx]; } + + void clear() { 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: + visual_weighting_array weights; + std::size_t count; +}; + 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 +111,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 +132,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 +158,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::ldexpf(1.0, -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 = 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,15 +237,16 @@ 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); } @@ -209,25 +255,26 @@ inline double csf_weight(double f, csf_model m, const csf_peak_t &pk) { // 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) { +inline visual_weights luma_visual_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). - 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}; + 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 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 + 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 - 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; + const float f_r = f_N * std::pow(2.0f, -(float)(lvl)) * std::sqrt(2.0f); + 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 @@ -237,60 +284,72 @@ inline std::vector luma_visual_weights(uint8_t dwt_levels, // --- 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; +inline 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 } - return std::vector(r, r + 15); + visual_weights result; + result.set_weights(r, 15); + return result; } -// 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 +// 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) { +inline visual_weights +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); } @@ -301,34 +360,35 @@ inline std::vector chroma_visual_weights(uint8_t dwt_levels, return luma_visual_weights(dwt_levels, vp); } - double sx = 1.0, sy = 1.0; // horizontal/vertical chroma subsampling factors + float sx = 1.0f, sy = 1.0f; // horizontal/vertical chroma subsampling factors if (chroma_format == 1) { // 4:2:0 - sx = 2.0; - sy = 2.0; + sx = 2.0f; + sy = 2.0f; } else if (chroma_format == 2) { // 4:2:2 - sx = 2.0; - sy = 1.0; + sx = 2.0f; + sy = 1.0f; } 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 + 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 - std::vector w; - w.reserve(static_cast(3) * dwt_levels); + visual_weights w; 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); + const float dx = f_Nx * std::pow(2.0f, -(float)(lvl)) * std::sqrt(2.0f); + const float dy = f_Ny * std::pow(2.0f, -(float)(lvl)) * std::sqrt(2.0f); + // 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 (vertical detail) - w.push_back(finite_weight(chroma_csf(dx, cp.a, cp.b))); // HL (horizontal detail) + 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; } From 1c1f9534d57e1e813f87f177ae6c1489b0ddf06d Mon Sep 17 00:00:00 2001 From: Aous Naman Date: Wed, 1 Jul 2026 20:11:16 +1000 Subject: [PATCH 2/8] Fixes compilation issues --- src/core/codestream/ojph_visual_weighting.h | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/core/codestream/ojph_visual_weighting.h b/src/core/codestream/ojph_visual_weighting.h index 2d288a7a..15aa283d 100644 --- a/src/core/codestream/ojph_visual_weighting.h +++ b/src/core/codestream/ojph_visual_weighting.h @@ -67,7 +67,6 @@ class visual_weights { std::size_t size() const { return count; } float operator[](std::size_t idx) { return weights[idx]; } - const float operator[](std::size_t idx) const { return weights[idx]; } void clear() { count = 0; } void push_back(float val) { @@ -184,7 +183,7 @@ inline q_scaling q_to_delta(uint8_t qfactor, uint8_t RI) { // 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 float eps0 = std::sqrt(0.5f) * std::ldexpf(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}; } @@ -228,7 +227,7 @@ inline csf_peak_t csf_peak(csf_model m) { auto scan = [](csf_model mm) { csf_peak_t pk{0.0f, 0.0f}; for (int i = 1; i < 100000; ++i) { - const float f = i * 0.002f; // 0.002 .. 200 cpd + 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; From 0c5873d07e63eeccfecc8e118d0f2a7798cd8381 Mon Sep 17 00:00:00 2001 From: Aous Naman Date: Thu, 2 Jul 2026 10:56:52 +1000 Subject: [PATCH 3/8] Removed dependence on std::array. --- src/core/codestream/ojph_visual_weighting.h | 24 +++++++++++---------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/src/core/codestream/ojph_visual_weighting.h b/src/core/codestream/ojph_visual_weighting.h index 15aa283d..bb20f131 100644 --- a/src/core/codestream/ojph_visual_weighting.h +++ b/src/core/codestream/ojph_visual_weighting.h @@ -43,30 +43,26 @@ // * 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 ///) +// 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 { -// Maximum number of DWT levels supported. -constexpr uint8_t max_dwt_levels = 32; -// Maximum number of visual weight entries (3 subbands per level). -constexpr size_t max_weight_size = 3 * max_dwt_levels; -// visual array type -using visual_weighting_array = std::array; // structure for our visual weights class visual_weights { public: - visual_weights() : weights(), count(0) {} + visual_weights() : weights{}, count(0) {} std::size_t size() const { return count; } - float operator[](std::size_t idx) { return weights[idx]; } + float operator[](std::size_t idx) { + assert(idx < max_weight_size); + return weights[idx]; + } void clear() { count = 0; } void push_back(float val) { @@ -82,8 +78,14 @@ class visual_weights { weights[i] = src[i]; } + public: + // 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: - visual_weighting_array weights; + float weights[max_weight_size]; std::size_t count; }; From 16f61c9f86e61433a4177e0ab8bdf27fa067c43e Mon Sep 17 00:00:00 2001 From: Aous Naman Date: Thu, 2 Jul 2026 11:08:19 +1000 Subject: [PATCH 4/8] Fixes a warning. --- src/core/codestream/ojph_params.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/core/codestream/ojph_params.cpp b/src/core/codestream/ojph_params.cpp index cc9d8c75..38671f82 100644 --- a/src/core/codestream/ojph_params.cpp +++ b/src/core/codestream/ojph_params.cpp @@ -1377,7 +1377,7 @@ namespace ojph { else { float w_b = (s == 0 || s > W_b.size()) ? - 1.0f : pow(W_b[W_b.size() - s], qfactor_power); + 1.0f : std::pow(W_b[W_b.size() - s], qfactor_power); delta_b = delta_ref / (w_g * w_b * G_c); } From bce9d5c7e679d815520dcf1ef65079822d470478 Mon Sep 17 00:00:00 2001 From: Aous Naman Date: Fri, 3 Jul 2026 13:07:14 +1000 Subject: [PATCH 5/8] Rearranged visual weighting and qfactor settings --- src/apps/ojph_compress/ojph_compress.cpp | 71 +---- src/core/codestream/ojph_params.cpp | 149 +++++++--- src/core/codestream/ojph_params_local.h | 18 +- src/core/codestream/ojph_visual_weighting.h | 306 ++++++++++---------- src/core/openjph/ojph_params.h | 134 +++++---- 5 files changed, 364 insertions(+), 314 deletions(-) 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 38671f82..abd2fe1c 100644 --- a/src/core/codestream/ojph_params.cpp +++ b/src/core/codestream/ojph_params.cpp @@ -418,6 +418,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) { @@ -1150,12 +1155,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) { @@ -1166,6 +1172,36 @@ 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); + comp_type t = ojph::param_qcd::ui8_2_comp_type(i < 3 ? i : 0); + 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 @@ -1192,7 +1228,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."); @@ -1222,7 +1259,8 @@ namespace ojph { } ////////////////////////////////////////////////////////////////////////// - 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."); @@ -1236,29 +1274,14 @@ namespace ojph { } ////////////////////////////////////////////////////////////////////////// - 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); - - param_qcd *p = get_qcc(comp_idx); - if (p == this) - p = add_qcc_object(comp_idx); + OJPH_ERROR(0x00050181, "Qfactor must be between 1 and 100, " + "but was set to %i.", qfactor); - p->qfactor = qfactor; - p->ctype = ctype; + this->qfactor = qfactor; } ////////////////////////////////////////////////////////////////////////// @@ -1317,55 +1340,48 @@ namespace ojph { // the following are used only when Qfactor is set float qfactor_power = 0; float delta_ref = 0; - float G_c = 1; + float G_c = 1.0f; const open_htj2k::visual_weighting_params vp; - open_htj2k::visual_weights W_b; - if (this->qfactor != QFACTOR_UNSET) { + open_htj2k::visual_weights v_weights; + open_htj2k::color_transform ct; + int chroma_factor = 0; + if (this->qfactor != QFACTOR_UNSET) + { const open_htj2k::q_scaling qs = open_htj2k::q_to_delta(this->qfactor, (ui8)ojph_min(16, 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; + ct = open_htj2k::resolve_color_transform(vp, this->is_color_trans); delta_ref = (float) (qs.delta_Q * open_htj2k::color_gain(ct, 0)); - G_c = (float) open_htj2k::color_gain(ct, comp_index); + G_c = (float) open_htj2k::color_gain(ct, this->ctype); // 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; + chroma_factor = 0; else if (this->sampling.x == 2 && this->sampling.y == 2) - sampling = 1; + chroma_factor = 1; else if (this->sampling.x == 2 && this->sampling.y == 1) - sampling = 2; + chroma_factor = 2; else 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 == 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); - } // 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; - + float w_g = 1.0f; + ui32 d = num_decomps - (ui32)(((int)s - 1) / 3); + ui32 sb = s != 0 ? (s - 1) % 3 + 1 : 0; if (num_decomps > 0) { //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 (s == 0) + if (sb == 0) { w_g = gain_l * gain_l; } - else if ((s - 1) % 3 == 2) + else if (sb == 3) { w_g = gain_h * gain_h; } else { w_g = gain_l * gain_h; } @@ -1376,8 +1392,18 @@ namespace ojph { delta_b = base_delta / w_g; else { - float w_b = (s == 0 || s > W_b.size()) ? - 1.0f : std::pow(W_b[W_b.size() - s], qfactor_power); + float w_b = 1.0f; + + bool result; + if (this->ctype == comp_type::OJPH_COMP_Y) + result = v_weights.luma_visual_weights(d - 1, sb, vp, w_b); + else + result = v_weights.chroma_visual_weights(d - 1, sb, vp, + this->ctype, chroma_factor, w_b, ct); + if (result == false) + OJPH_ERROR(0x00050162, + "Something is wrong with visual weight settings"); + w_b = std::pow(w_b, qfactor_power); delta_b = delta_ref / (w_g * w_b * G_c); } @@ -1788,6 +1814,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 4651a4ab..570e037d 100644 --- a/src/core/codestream/ojph_params_local.h +++ b/src/core/codestream/ojph_params_local.h @@ -696,6 +696,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); } @@ -712,11 +715,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, @@ -730,6 +734,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); @@ -747,7 +753,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; @@ -798,7 +804,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 bb20f131..560c17ea 100644 --- a/src/core/codestream/ojph_visual_weighting.h +++ b/src/core/codestream/ojph_visual_weighting.h @@ -53,41 +53,6 @@ namespace open_htj2k { -// structure for our visual weights -class visual_weights { - public: - visual_weights() : weights{}, count(0) {} - - std::size_t size() const { return count; } - float operator[](std::size_t idx) { - assert(idx < max_weight_size); - return weights[idx]; - } - - void clear() { 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]; - } - - public: - // 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; -}; enum class csf_model { legacy_table, // exact Zeng Table 2 weights (default; bit-identical output) @@ -251,38 +216,6 @@ inline float csf_weight(float f, csf_model m, const csf_peak_t &pk) { 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 visual_weights luma_visual_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)) * std::sqrt(2.0f); - 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; -} - // --- chroma (chrominance) analytic weighting ------------------------------- // // Chrominance CSF is low-pass (not band-pass like luminance) and rolls off @@ -310,88 +243,171 @@ inline chroma_csf_params chroma_params_for(int comp_index) { 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 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 +// structure for our visual weights +class visual_weights { + public: + visual_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. + bool luma_visual_weights(uint8_t dwt_level, uint8_t orientation, + const visual_weighting_params &vp, float &weight) + { + if (vp.model == csf_model::legacy_table) + return legacy_luma_lookup(dwt_level, orientation, weight); + else + return luma_lookup(dwt_level, orientation, vp, weight); } - visual_weights result; - result.set_weights(r, 15); - return result; -} -// 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 visual_weights -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); + // 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. + bool chroma_visual_weights(uint8_t dwt_level, uint8_t orientation, + const visual_weighting_params &vp, int comp_index, int chroma_format, + float &weight, color_transform ct = color_transform::ict) + { + if (vp.model == csf_model::legacy_table) + return legacy_chroma_lookup(comp_index, chroma_format, dwt_level, + orientation, weight); + + if (ct == color_transform::none) + return legacy_luma_lookup(dwt_level, orientation, weight); + else + return chroma_lookup(comp_index, chroma_format, dwt_level, + orientation, vp, weight); } - float sx = 1.0f, sy = 1.0f; // horizontal/vertical 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; - } + private: + //////////////////////////////////////// legacy tables + bool legacy_luma_lookup(uint8_t dwt_level, uint8_t orientation, + float &weight) + { + return legacy_look_up(dwt_level, orientation, luma, weight); + } - 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)) * std::sqrt(2.0f); - const float dy = f_Ny * std::pow(2.0f, -(float)(lvl)) * std::sqrt(2.0f); - // 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; -} + bool legacy_chroma_lookup(int comp_index, int chroma_format, + uint8_t dwt_level, uint8_t orientation, float &weight) + { + 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 + } + return legacy_look_up(dwt_level, orientation, r, weight); + } + + bool legacy_look_up(uint8_t dwt_level, uint8_t orientation, + const float *tbl, float &weight) + { + if (orientation > 3) + { weight = 0.0f; return false; } + if (orientation == 0 || dwt_level >= 5) + { weight = 1.0f; return true; } + else + { weight = tbl[dwt_level * 3 + 3 - orientation]; return true; } + } + + //////////////////////////////////////// formula based + bool luma_lookup(uint8_t dwt_level, uint8_t orientation, + const visual_weighting_params &vp, float &weight) + { + if (orientation >= 3) { weight = 0.0f; return false; } + if (orientation == 0) { weight = 1.0f; return true; } + + 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 + + // Geometric-mean radial center of octave band + // [f_N/2^dwt_level, f_N/2^(dwt_level-1)]. + const float f_r = + f_N * std::pow(2.0f, -(float)(dwt_level)) * (float)M_SQRT2; + const float f_hh = f_r * vp.hh_factor; + if (orientation == 1) + weight = csf_weight(f_r, vp.model, pk); + else if (orientation == 2) + weight = csf_weight(f_r, vp.model, pk); + else if (orientation == 3) + weight = csf_weight(f_hh, vp.model, pk); + + return true; + } + + bool chroma_lookup(int comp_index, int chroma_format, uint8_t dwt_level, + uint8_t orientation, const visual_weighting_params &vp, float &weight) + { + if (orientation >= 3) { weight = 0.0f; return false; } + if (orientation == 0) { weight = 1.0f; return true; } + + float sx = 1.0f, sy = 1.0f; // horz/vert chroma subsampling factors + if (chroma_format == 1) + sx = sy = 2.0f; // 4:2:0 + 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 + + const float t = std::pow(2.0f, -(float)(dwt_level)) * (float)M_SQRT2; + const float dx = f_Nx * t; + const float dy = f_Ny * t; + + // 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. + + if (orientation == 1) + weight = finite_weight(chroma_csf(dx, cp.a, cp.b)); + else if (orientation == 2) + weight = finite_weight(chroma_csf(dy, cp.a, cp.b)); + else if (orientation == 3) + weight = finite_weight(chroma_csf(f_hh, cp.a, cp.b)); + + return true; + } + + 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; + // Luminance and Chromium legacy weights + const float luma[15] = + {0.0901f, 0.2758f, 0.2758f, 0.7018f, 0.8378f, 0.8378f, 1.0000f, 1.0000f, + 1.0000f, 1.0000f, 1.0000f, 1.0000f, 1.0000f, 1.0000f, 1.0000f}; + 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}; + 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}; + 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.0000f, 1.0000f, 1.0000f}; + 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.0000f, 1.0000f, 1.0000f}; + 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.0000f}; + 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.0000f}; +}; } // namespace open_htj2k diff --git a/src/core/openjph/ojph_params.h b/src/core/openjph/ojph_params.h index dce330d2..785289ce 100644 --- a/src/core/openjph/ojph_params.h +++ b/src/core/openjph/ojph_params.h @@ -2,21 +2,21 @@ // This software is released under the 2-Clause BSD license, included // below. // -// Copyright (c) 2019, Aous Naman +// 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 @@ -157,50 +157,74 @@ namespace ojph { /***************************************************************************/ /** * @brief Quantization parameters object - * + * */ 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 ui8 comp_type_2_ui8(comp_type ctype) + { return static_cast(ctype); } + 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) {} /** - * @brief Set the irreversible quantization base delta. - * - * This represents the default base delta and influences QCD marker + * @brief Set the irreversible quantization base delta. + * + * This represents the default base delta and influences QCD marker * segment - * - * @param delta + * + * @param delta */ void set_irrev_quant(float delta); /** - * @brief Set the irreversible quantization base delta for a specific + * @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 - * - * This represents the default base delta for component comp_idx, and + * + * This represents the default base delta for component comp_idx, and * influences QCC marker segment for the component, inserting one * if needed, which is usually the case. - * - * @param comp_idx - * @param delta + * + * @param comp_idx + * @param delta */ void set_irrev_quant(ui32 comp_idx, float delta); /** * @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); @@ -212,82 +236,82 @@ namespace ojph { /** * @brief non-linearity point transformation object * (implements NLT marker segment) - * - * There are a few things to know here. - * The NLT marker segment contains the nonlinearity type and the + * + * There are a few things to know here. + * The NLT marker segment contains the nonlinearity type and the * bit depth and signedness of the component to which it applies. - * There is the default component ALL_COMPS which applies to all + * There is the default component ALL_COMPS which applies to all * components unless it is overridden by another NLT segment marker. * The library checks that the settings make sense, and also make * sure that bit depth and signedness are correct, creating any missing * NLT marker segments in the process. * If all components have the same bit depth and signedness, and need - * nonlinearity type 3 (Binary Complement to Sign Magnitude Conversion), + * nonlinearity type 3 (Binary Complement to Sign Magnitude Conversion), * then the best option is to set ALL_COMPS to type 3. - * Otherwise, the best option is to set type 3 only to components that + * Otherwise, the best option is to set type 3 only to components that * need it, leaving out the default ALL_COMPS nonlinearity not set. - * Another option is for the end-user can set the ALL_COMPS to type 3, - * and then put exception for the components that does not need type 3, + * Another option is for the end-user can set the ALL_COMPS to type 3, + * and then put exception for the components that does not need type 3, * by setting them to type 0. - * + * * The library, during validity check, which is run when the codestream * is created for writing, will do the following: - * -- If ALL_COMPS is set to type 0, it will be ignored, and the + * -- If ALL_COMPS is set to type 0, it will be ignored, and the * codestream will NOT have the corresponding NLT marker segment. * -- If ALL_COMPS is set to type 3, then the following will happen: - * - If all the components (except those with type 0 set for them) have - * the same bit depth and signedness, then the ALL_COMPS NLT marker + * - If all the components (except those with type 0 set for them) have + * the same bit depth and signedness, then the ALL_COMPS NLT marker * segment will be respected and inserted into the codestream. * Of course, components with NLT 0 will also have the corresponding * NLT marker segment inserted. * - If components, for which no NTL type 0 is specified, have differing - * bit depth or signedness, then the ALL_COMPS will be ignored, and + * bit depth or signedness, then the ALL_COMPS will be ignored, and * NLT markers are inserted for each component that needs type 3. * Components that have their component field larger than the number of * components in the codestream are removed. - * - * It also worth noting that type 3 nonlinearity has no effect on - * positive image samples. It is also not recommended for integer-valued - * types. It is only recommended for floating-point image samples, for - * which some of the samples are negative, where type 3 nonlinearity - * should be beneficial. This is because the encoding engine expects - * two-complement representation for negative values while floating point - * numbers have a sign bit followed by an exponent, which has a biased + * + * It also worth noting that type 3 nonlinearity has no effect on + * positive image samples. It is also not recommended for integer-valued + * types. It is only recommended for floating-point image samples, for + * which some of the samples are negative, where type 3 nonlinearity + * should be beneficial. This is because the encoding engine expects + * two-complement representation for negative values while floating point + * numbers have a sign bit followed by an exponent, which has a biased * integer representation. The core idea is to make floating-point * representation more compatible with integer representation. - * + * */ class OJPH_EXPORT param_nlt { public: enum special_comp_num : ui16 { ALL_COMPS = 65535 }; - enum nonlinearity : ui8 { + enum nonlinearity : ui8 { OJPH_NLT_NO_NLT = 0, // supported OJPH_NLT_GAMMA_STYLE_NLT = 1, // not supported OJPH_NLT_LUT_STYLE_NLT = 2, // not supported OJPH_NLT_BINARY_COMPLEMENT_NLT = 3, // supported - OJPH_NLT_UNDEFINED = 255 // This is used internally and is - // not part of the standard + OJPH_NLT_UNDEFINED = 255 // This is used internally and is + // not part of the standard }; public: param_nlt(local::param_nlt* p) : state(p) {} /** - * @brief enables or disables type 3 nonlinearity for a component + * @brief enables or disables type 3 nonlinearity for a component * or the default setting - * + * * When creating a codestream for writing, call this function before * you call codestream::write_headers. - * - * + * + * * @param comp_num: component number, or 65535 for the default setting * @param type: desired non-linearity from enum nonlinearity */ void set_nonlinear_transform(ui32 comp_num, ui8 nl_type); /** - * @brief get the nonlinearity type associated with comp_num, which + * @brief get the nonlinearity type associated with comp_num, which * should be one from enum nonlinearity * * @param comp_num: component number, or 65535 for the default setting @@ -296,7 +320,7 @@ namespace ojph { * @param type: nonlinearity type * @return true if the nonlinearity for comp_num is set */ - bool get_nonlinear_transform(ui32 comp_num, ui8& bit_depth, + bool get_nonlinear_transform(ui32 comp_num, ui8& bit_depth, bool& is_signed, ui8& nl_type) const; private: From 05673eae4f591ec78ad98a9ed457df7e40f8644d Mon Sep 17 00:00:00 2001 From: Aous Naman Date: Fri, 3 Jul 2026 13:18:33 +1000 Subject: [PATCH 6/8] Fixes warnings. --- src/core/codestream/ojph_params.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/core/codestream/ojph_params.cpp b/src/core/codestream/ojph_params.cpp index 50fa37d9..e6e7433d 100644 --- a/src/core/codestream/ojph_params.cpp +++ b/src/core/codestream/ojph_params.cpp @@ -1201,7 +1201,8 @@ namespace ojph { if (q == this) { q = add_qcc_object(i); - comp_type t = ojph::param_qcd::ui8_2_comp_type(i < 3 ? i : 0); + ui8 ci = (ui8)(i < 3u ? i : 0u); + comp_type t = ojph::param_qcd::ui8_2_comp_type(ci); set_qfactor(i, t, this->qfactor); } } @@ -1377,8 +1378,8 @@ namespace ojph { { // compute square root of the enery gain factor W_g float w_g = 1.0f; - ui32 d = num_decomps - (ui32)(((int)s - 1) / 3); - ui32 sb = s != 0 ? (s - 1) % 3 + 1 : 0; + ui8 d = (ui8)(num_decomps - (ui32)(((int)s - 1) / 3)); + ui8 sb = (ui8)(s != 0 ? (s - 1) % 3 + 1 : 0); if (num_decomps > 0) { //In C++, division result truncates towards zero From 8b543b0ac64d42672faa6ba9ff94f75fe83ed56c Mon Sep 17 00:00:00 2001 From: Aous Naman Date: Fri, 3 Jul 2026 17:51:47 +1000 Subject: [PATCH 7/8] Reverts visual weights to an older version and create a factory function for it. Also break setting SPcod into three functions. --- src/core/codestream/ojph_params.cpp | 147 +++++----- src/core/codestream/ojph_params_local.h | 2 + src/core/codestream/ojph_visual_weighting.h | 282 ++++++++++---------- src/core/openjph/ojph_params.h | 2 - 4 files changed, 231 insertions(+), 202 deletions(-) diff --git a/src/core/codestream/ojph_params.cpp b/src/core/codestream/ojph_params.cpp index e6e7433d..f09f0f53 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 { //////////////////////////////////////////////////////////////////////////// @@ -1261,7 +1262,10 @@ 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); } } @@ -1277,7 +1281,6 @@ 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(); - } ////////////////////////////////////////////////////////////////////////// @@ -1344,87 +1347,111 @@ 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.0f; - const open_htj2k::visual_weighting_params vp; - open_htj2k::visual_weights v_weights; - open_htj2k::color_transform ct; - int chroma_factor = 0; - if (this->qfactor != QFACTOR_UNSET) - { - const open_htj2k::q_scaling qs = open_htj2k::q_to_delta(this->qfactor, - (ui8)ojph_min(16, this->bit_depth)); - qfactor_power = (float) qs.qfactor_power; - ct = open_htj2k::resolve_color_transform(vp, this->is_color_trans); - delta_ref = (float) (qs.delta_Q * open_htj2k::color_gain(ct, 0)); - G_c = (float) open_htj2k::color_gain(ct, this->ctype); - - // chroma_format 0 = 4:4:4, 1 = 4:2:0, 2 = 4:2:2 - if (this->sampling.x == 1 && this->sampling.y == 1) - chroma_factor = 0; - else if (this->sampling.x == 2 && this->sampling.y == 2) - chroma_factor = 1; - else if (this->sampling.x == 2 && this->sampling.y == 1) - chroma_factor = 2; - else - OJPH_ERROR(0x00050161, "Qfactor can only be used on components " - "with 4:4:4, 4:2:2 or 4:2:0 sampling"); - } - // 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.0f; - ui8 d = (ui8)(num_decomps - (ui32)(((int)s - 1) / 3)); - ui8 sb = (ui8)(s != 0 ? (s - 1) % 3 + 1 : 0); + float w_g = 1.0; + if (num_decomps > 0) { //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 (sb == 0) + if (s == 0) { w_g = gain_l * gain_l; } - else if (sb == 3) + else if ((s - 1) % 3 == 2) { w_g = gain_h * gain_h; } else { w_g = gain_l * gain_h; } } - float delta_b; - if (this->qfactor == QFACTOR_UNSET) - delta_b = base_delta / w_g; - else + 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++) + { + // compute square root of the enery gain factor W_g + float w_g = 1.0; + + if (num_decomps > 0) { - float w_b = 1.0f; + //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); - bool result; - if (this->ctype == comp_type::OJPH_COMP_Y) - result = v_weights.luma_visual_weights(d - 1, sb, vp, w_b); + if (s == 0) + { w_g = gain_l * gain_l; } + else if ((s - 1) % 3 == 2) + { w_g = gain_h * gain_h; } else - result = v_weights.chroma_visual_weights(d - 1, sb, vp, - this->ctype, chroma_factor, w_b, ct); - if (result == false) - OJPH_ERROR(0x00050162, - "Something is wrong with visual weight settings"); - w_b = std::pow(w_b, qfactor_power); - delta_b = delta_ref / (w_g * w_b * G_c); + { w_g = gain_l * gain_h; } } - 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 delta_b; + 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); + printf("%f\n", w_b); + 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 { diff --git a/src/core/codestream/ojph_params_local.h b/src/core/codestream/ojph_params_local.h index cb172c47..ab25fa8b 100644 --- a/src/core/codestream/ojph_params_local.h +++ b/src/core/codestream/ojph_params_local.h @@ -786,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); diff --git a/src/core/codestream/ojph_visual_weighting.h b/src/core/codestream/ojph_visual_weighting.h index 560c17ea..33c7870c 100644 --- a/src/core/codestream/ojph_visual_weighting.h +++ b/src/core/codestream/ojph_visual_weighting.h @@ -27,7 +27,8 @@ // 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. @@ -53,7 +54,6 @@ namespace open_htj2k { - enum class csf_model { legacy_table, // exact Zeng Table 2 weights (default; bit-identical output) mannos_sakrison, // analytic Mannos-Sakrison CSF @@ -245,140 +245,158 @@ inline chroma_csf_params chroma_params_for(int comp_index) { // structure for our visual weights class visual_weights { - public: - visual_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. - bool luma_visual_weights(uint8_t dwt_level, uint8_t orientation, - const visual_weighting_params &vp, float &weight) - { - if (vp.model == csf_model::legacy_table) - return legacy_luma_lookup(dwt_level, orientation, weight); - else - return luma_lookup(dwt_level, orientation, vp, weight); - } - - // 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. - bool chroma_visual_weights(uint8_t dwt_level, uint8_t orientation, - const visual_weighting_params &vp, int comp_index, int chroma_format, - float &weight, color_transform ct = color_transform::ict) - { - if (vp.model == csf_model::legacy_table) - return legacy_chroma_lookup(comp_index, chroma_format, dwt_level, - orientation, weight); - - if (ct == color_transform::none) - return legacy_luma_lookup(dwt_level, orientation, weight); - else - return chroma_lookup(comp_index, chroma_format, dwt_level, - orientation, vp, weight); - } - - private: - //////////////////////////////////////// legacy tables - bool legacy_luma_lookup(uint8_t dwt_level, uint8_t orientation, - float &weight) - { - return legacy_look_up(dwt_level, orientation, luma, weight); - } - bool legacy_chroma_lookup(int comp_index, int chroma_format, - uint8_t dwt_level, uint8_t orientation, float &weight) - { - 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 - } - return legacy_look_up(dwt_level, orientation, r, weight); - } - - bool legacy_look_up(uint8_t dwt_level, uint8_t orientation, - const float *tbl, float &weight) + 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 (orientation > 3) - { weight = 0.0f; return false; } - if (orientation == 0 || dwt_level >= 5) - { weight = 1.0f; return true; } + if (comp_index == 0) + return make_luma_weights(dwt_levels, vp); else - { weight = tbl[dwt_level * 3 + 3 - orientation]; return true; } + return + make_chroma_weights(dwt_levels, vp, comp_index, chroma_format, ct); } - //////////////////////////////////////// formula based - bool luma_lookup(uint8_t dwt_level, uint8_t orientation, - const visual_weighting_params &vp, float &weight) + const float* get_weights() const { return weights; } + size_t get_count() const { return count; } + + 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) { - if (orientation >= 3) { weight = 0.0f; return false; } - if (orientation == 0) { weight = 1.0f; return true; } + 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 - // Geometric-mean radial center of octave band - // [f_N/2^dwt_level, f_N/2^(dwt_level-1)]. - const float f_r = - f_N * std::pow(2.0f, -(float)(dwt_level)) * (float)M_SQRT2; - const float f_hh = f_r * vp.hh_factor; - if (orientation == 1) - weight = csf_weight(f_r, vp.model, pk); - else if (orientation == 2) - weight = csf_weight(f_r, vp.model, pk); - else if (orientation == 3) - weight = csf_weight(f_hh, vp.model, pk); - - return true; + 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; } - bool chroma_lookup(int comp_index, int chroma_format, uint8_t dwt_level, - uint8_t orientation, const visual_weighting_params &vp, float &weight) + // 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 (orientation >= 3) { weight = 0.0f; return false; } - if (orientation == 0) { weight = 1.0f; return true; } + 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) - sx = sy = 2.0f; // 4:2:0 - else if (chroma_format == 2) // 4:2:2 - { sx = 2.0f; sy = 1.0f; } + 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 - - const float t = std::pow(2.0f, -(float)(dwt_level)) * (float)M_SQRT2; - const float dx = f_Nx * t; - const float dy = f_Ny * t; - - // 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. - - if (orientation == 1) - weight = finite_weight(chroma_csf(dx, cp.a, cp.b)); - else if (orientation == 2) - weight = finite_weight(chroma_csf(dy, cp.a, cp.b)); - else if (orientation == 3) - weight = finite_weight(chroma_csf(f_hh, cp.a, cp.b)); - - return true; + 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: @@ -386,28 +404,12 @@ class visual_weights { 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; - // Luminance and Chromium legacy weights - const float luma[15] = - {0.0901f, 0.2758f, 0.2758f, 0.7018f, 0.8378f, 0.8378f, 1.0000f, 1.0000f, - 1.0000f, 1.0000f, 1.0000f, 1.0000f, 1.0000f, 1.0000f, 1.0000f}; - 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}; - 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}; - 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.0000f, 1.0000f, 1.0000f}; - 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.0000f, 1.0000f, 1.0000f}; - 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.0000f}; - 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.0000f}; + + 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 fecfa8fa..4d0f284b 100644 --- a/src/core/openjph/ojph_params.h +++ b/src/core/openjph/ojph_params.h @@ -175,8 +175,6 @@ namespace ojph { OJPH_COMP_CR = 2, OJPH_COMP_UNDEFINED = 0xFF }; - static ui8 comp_type_2_ui8(comp_type ctype) - { return static_cast(ctype); } static comp_type ui8_2_comp_type(ui8 c) { if (c >= OJPH_COMP_Y && c <= OJPH_COMP_CR) From 771eb0ef4fa7d000570835fc56fcbfa7855fcb7e Mon Sep 17 00:00:00 2001 From: Aous Naman Date: Fri, 3 Jul 2026 17:57:41 +1000 Subject: [PATCH 8/8] Removes one forgotten line --- src/core/codestream/ojph_params.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/core/codestream/ojph_params.cpp b/src/core/codestream/ojph_params.cpp index f09f0f53..b87106df 100644 --- a/src/core/codestream/ojph_params.cpp +++ b/src/core/codestream/ojph_params.cpp @@ -1434,7 +1434,6 @@ namespace ojph { if (s > 0 && s <= count) w_b = std::pow(w[count - s], qfactor_power); delta_b = delta_ref / (w_g * w_b * G_c); - printf("%f\n", w_b); encode_SPqcd(s, delta_b); } }