From d071eca040ec317dded77e6c741ccd27db5db799 Mon Sep 17 00:00:00 2001 From: Sam Li Date: Tue, 17 Mar 2026 20:23:09 -0700 Subject: [PATCH 1/3] move function msb_position() to sperr_helper --- include/SPECK3D_INT_ENC.h | 2 - include/sperr_helper.h | 4 ++ src/SPECK3D_INT_ENC.cpp | 69 +++++++++---------------- src/sperr_helper.cpp | 23 +++++++++ test_scripts/sperr_helper_unit_test.cpp | 31 +++++++++++ 5 files changed, 83 insertions(+), 46 deletions(-) diff --git a/include/SPECK3D_INT_ENC.h b/include/SPECK3D_INT_ENC.h index c078147b..f1faff67 100644 --- a/include/SPECK3D_INT_ENC.h +++ b/include/SPECK3D_INT_ENC.h @@ -47,8 +47,6 @@ class SPECK3D_INT_ENC final : public SPECK3D_INT { // m_bitplane_init(). Significance tests compare m_morton_buf entries against this value. int8_t m_morton_threshold = -1; void m_deposit_set(Set3D); - // Returns the bit position of the most significant bit (0-based), or -1 for zero. - auto m_msb_position(uint_type v) const -> int8_t; }; }; // namespace sperr diff --git a/include/sperr_helper.h b/include/sperr_helper.h index c9385f0c..7e90a58a 100644 --- a/include/sperr_helper.h +++ b/include/sperr_helper.h @@ -166,6 +166,10 @@ auto calc_stats(const T* arr1, const T* arr2, size_t arr_len, size_t omp_nthread template auto kahan_summation(const T*, size_t) -> T; +// Returns the bit position of the most significant bit (0-based), or -1 for zero. +template +auto msb_position(T v) -> int8_t; + // Given a whole volume size and a desired chunk size, this helper function // returns a list of chunks specified by 6 integers: // chunk[0], [2], [4]: starting index of this chunk in X, Y, and Z; diff --git a/src/SPECK3D_INT_ENC.cpp b/src/SPECK3D_INT_ENC.cpp index b698b4a3..90336bcf 100644 --- a/src/SPECK3D_INT_ENC.cpp +++ b/src/SPECK3D_INT_ENC.cpp @@ -5,10 +5,6 @@ #include // std::memcpy() #include -#if __cplusplus >= 202002L -#include -#endif - template void sperr::SPECK3D_INT_ENC::m_deposit_set(Set3D set) { @@ -17,7 +13,7 @@ void sperr::SPECK3D_INT_ENC::m_deposit_set(Set3D set) return; case 1: { auto id = set.start_z * m_dims[0] * m_dims[1] + set.start_y * m_dims[0] + set.start_x; - m_morton_buf[set.morton_idx] = m_msb_position(m_coeff_buf[id]); + m_morton_buf[set.morton_idx] = sperr::msb_position(m_coeff_buf[id]); return; } case 2: { @@ -26,7 +22,7 @@ void sperr::SPECK3D_INT_ENC::m_deposit_set(Set3D set) // Deposit the 1st element. auto id = set.start_z * m_dims[0] * m_dims[1] + set.start_y * m_dims[0] + set.start_x; auto morton_id = set.morton_idx; - m_morton_buf[morton_id] = m_msb_position(m_coeff_buf[id]); + m_morton_buf[morton_id] = sperr::msb_position(m_coeff_buf[id]); // Deposit the 2nd element. if (set.length_x == 2) @@ -35,7 +31,7 @@ void sperr::SPECK3D_INT_ENC::m_deposit_set(Set3D set) id += m_dims[0]; else id += m_dims[0] * m_dims[1]; - m_morton_buf[++morton_id] = m_msb_position(m_coeff_buf[id]); + m_morton_buf[++morton_id] = sperr::msb_position(m_coeff_buf[id]); return; } @@ -45,51 +41,51 @@ void sperr::SPECK3D_INT_ENC::m_deposit_set(Set3D set) if (set.length_x == 2 && set.length_y == 2) { // Element (0, 0, 0) - m_morton_buf[morton_id] = m_msb_position(m_coeff_buf[id]); + m_morton_buf[morton_id] = sperr::msb_position(m_coeff_buf[id]); // Element (1, 0, 0) - m_morton_buf[++morton_id] = m_msb_position(m_coeff_buf[id + 1]); + m_morton_buf[++morton_id] = sperr::msb_position(m_coeff_buf[id + 1]); // Element (0, 1, 0) auto id2 = id + m_dims[0]; - m_morton_buf[++morton_id] = m_msb_position(m_coeff_buf[id2]); + m_morton_buf[++morton_id] = sperr::msb_position(m_coeff_buf[id2]); // Element (1, 1, 0) - m_morton_buf[++morton_id] = m_msb_position(m_coeff_buf[++id2]); + m_morton_buf[++morton_id] = sperr::msb_position(m_coeff_buf[++id2]); return; } else if (set.length_x == 2 && set.length_z == 2) { // Element (0, 0, 0) - m_morton_buf[morton_id] = m_msb_position(m_coeff_buf[id]); + m_morton_buf[morton_id] = sperr::msb_position(m_coeff_buf[id]); // Element (1, 0, 0) - m_morton_buf[++morton_id] = m_msb_position(m_coeff_buf[id + 1]); + m_morton_buf[++morton_id] = sperr::msb_position(m_coeff_buf[id + 1]); // Element (0, 0, 1) auto id2 = id + m_dims[0] * m_dims[1]; - m_morton_buf[++morton_id] = m_msb_position(m_coeff_buf[id2]); + m_morton_buf[++morton_id] = sperr::msb_position(m_coeff_buf[id2]); // Element (1, 0, 1) - m_morton_buf[++morton_id] = m_msb_position(m_coeff_buf[++id2]); + m_morton_buf[++morton_id] = sperr::msb_position(m_coeff_buf[++id2]); return; } else if (set.length_y == 2 && set.length_z == 2) { // Element (0, 0, 0) - m_morton_buf[morton_id] = m_msb_position(m_coeff_buf[id]); + m_morton_buf[morton_id] = sperr::msb_position(m_coeff_buf[id]); // Element (0, 1, 0) auto id2 = id + m_dims[0]; - m_morton_buf[++morton_id] = m_msb_position(m_coeff_buf[id2]); + m_morton_buf[++morton_id] = sperr::msb_position(m_coeff_buf[id2]); // Element (0, 0, 1) id2 = id + m_dims[0] * m_dims[1]; - m_morton_buf[++morton_id] = m_msb_position(m_coeff_buf[id2]); + m_morton_buf[++morton_id] = sperr::msb_position(m_coeff_buf[id2]); // Element (0, 1, 1) id2 += m_dims[0]; - m_morton_buf[++morton_id] = m_msb_position(m_coeff_buf[id2]); + m_morton_buf[++morton_id] = sperr::msb_position(m_coeff_buf[id2]); return; } @@ -101,31 +97,31 @@ void sperr::SPECK3D_INT_ENC::m_deposit_set(Set3D set) // Element (0, 0, 0) const auto id = set.start_z * m_dims[0] * m_dims[1] + set.start_y * m_dims[0] + set.start_x; auto morton_id = set.morton_idx; - m_morton_buf[morton_id] = m_msb_position(m_coeff_buf[id]); + m_morton_buf[morton_id] = sperr::msb_position(m_coeff_buf[id]); // Element (1, 0, 0) - m_morton_buf[++morton_id] = m_msb_position(m_coeff_buf[id + 1]); + m_morton_buf[++morton_id] = sperr::msb_position(m_coeff_buf[id + 1]); // Element (0, 1, 0) auto id2 = id + m_dims[0]; - m_morton_buf[++morton_id] = m_msb_position(m_coeff_buf[id2]); + m_morton_buf[++morton_id] = sperr::msb_position(m_coeff_buf[id2]); // Element (1, 1, 0) - m_morton_buf[++morton_id] = m_msb_position(m_coeff_buf[++id2]); + m_morton_buf[++morton_id] = sperr::msb_position(m_coeff_buf[++id2]); // Element (0, 0, 1) id2 = id + m_dims[0] * m_dims[1]; - m_morton_buf[++morton_id] = m_msb_position(m_coeff_buf[id2]); + m_morton_buf[++morton_id] = sperr::msb_position(m_coeff_buf[id2]); // Element (1, 0, 1) - m_morton_buf[++morton_id] = m_msb_position(m_coeff_buf[++id2]); + m_morton_buf[++morton_id] = sperr::msb_position(m_coeff_buf[++id2]); // Element (0, 1, 1) id2 = id + m_dims[0] * (m_dims[1] + 1); - m_morton_buf[++morton_id] = m_msb_position(m_coeff_buf[id2]); + m_morton_buf[++morton_id] = sperr::msb_position(m_coeff_buf[id2]); // Element (1, 1, 1) - m_morton_buf[++morton_id] = m_msb_position(m_coeff_buf[++id2]); + m_morton_buf[++morton_id] = sperr::msb_position(m_coeff_buf[++id2]); return; } @@ -189,7 +185,7 @@ void sperr::SPECK3D_INT_ENC::m_process_P(size_t idx, size_t morton, size_t& c bool is_sig = true; if (output) { - assert(m_msb_position(m_coeff_buf[idx]) == m_morton_buf[morton]); + assert(sperr::msb_position(m_coeff_buf[idx]) == m_morton_buf[morton]); is_sig = (m_morton_buf[morton] >= m_morton_threshold); m_bit_buffer.wbit(is_sig); } @@ -215,25 +211,10 @@ void sperr::SPECK3D_INT_ENC::m_process_P_lite(size_t idx) } } -template -auto sperr::SPECK3D_INT_ENC::m_msb_position(uint_type v) const -> int8_t -{ -#if __cplusplus >= 202002L - return static_cast(sizeof(uint_type) * 8 - 1 - std::countl_zero(v)); -#else - int8_t pos = -1; - while (v) { - v >>= 1; - pos++; - } - return pos; -#endif -} - template void sperr::SPECK3D_INT_ENC::m_bitplane_init() { - m_morton_threshold = m_msb_position(m_threshold); + m_morton_threshold = sperr::msb_position(m_threshold); } template diff --git a/src/sperr_helper.cpp b/src/sperr_helper.cpp index a2543834..b72aae57 100644 --- a/src/sperr_helper.cpp +++ b/src/sperr_helper.cpp @@ -7,6 +7,10 @@ #include #include +#if __cplusplus >= 202002L +#include +#endif + #ifdef __AVX2__ #include #endif @@ -641,3 +645,22 @@ auto sperr::calc_mean_var(const T* arr, size_t len, size_t omp_nthreads) -> std: } template auto sperr::calc_mean_var(const float*, size_t, size_t) -> std::array; template auto sperr::calc_mean_var(const double*, size_t, size_t) -> std::array; + +template +auto sperr::msb_position(T v) -> int8_t +{ +#if __cplusplus >= 202002L + return static_cast(sizeof(T) * 8 - 1 - std::countl_zero(v)); +#else + int8_t pos = -1; + while (v) { + v >>= 1; + pos++; + } + return pos; +#endif +} +template auto sperr::msb_position(uint8_t) -> int8_t; +template auto sperr::msb_position(uint16_t) -> int8_t; +template auto sperr::msb_position(uint32_t) -> int8_t; +template auto sperr::msb_position(uint64_t) -> int8_t; diff --git a/test_scripts/sperr_helper_unit_test.cpp b/test_scripts/sperr_helper_unit_test.cpp index ef273be0..b8bccead 100644 --- a/test_scripts/sperr_helper_unit_test.cpp +++ b/test_scripts/sperr_helper_unit_test.cpp @@ -293,4 +293,35 @@ TEST(sperr_helper, read_sections) EXPECT_EQ(buf, buf2); } +TEST(sperr_helper, msb_position) +{ + // Zero returns -1 + EXPECT_EQ(sperr::msb_position(uint8_t{0}), -1); + EXPECT_EQ(sperr::msb_position(uint16_t{0}), -1); + EXPECT_EQ(sperr::msb_position(uint32_t{0}), -1); + EXPECT_EQ(sperr::msb_position(uint64_t{0}), -1); + + // Powers of 2 + EXPECT_EQ(sperr::msb_position(uint32_t{1}), 0); + EXPECT_EQ(sperr::msb_position(uint32_t{2}), 1); + EXPECT_EQ(sperr::msb_position(uint32_t{4}), 2); + EXPECT_EQ(sperr::msb_position(uint32_t{1024}), 10); + + // Non-powers of 2 + EXPECT_EQ(sperr::msb_position(uint32_t{3}), 1); + EXPECT_EQ(sperr::msb_position(uint32_t{5}), 2); + EXPECT_EQ(sperr::msb_position(uint32_t{255}), 7); + EXPECT_EQ(sperr::msb_position(uint32_t{256}), 8); + + // Max values for each type + EXPECT_EQ(sperr::msb_position(std::numeric_limits::max()), 7); + EXPECT_EQ(sperr::msb_position(std::numeric_limits::max()), 15); + EXPECT_EQ(sperr::msb_position(std::numeric_limits::max()), 31); + EXPECT_EQ(sperr::msb_position(std::numeric_limits::max()), 63); + + // Large 64-bit values + EXPECT_EQ(sperr::msb_position(uint64_t{1} << 40), 40); + EXPECT_EQ(sperr::msb_position(uint64_t{1} << 62), 62); +} + } // namespace From 6c5e33ef9a3eb86f6b3463e710990f5bcd1f538e Mon Sep 17 00:00:00 2001 From: Sam Li Date: Tue, 17 Mar 2026 20:59:08 -0700 Subject: [PATCH 2/3] SPECK2D_INT_ENC also uses an int8_t buffer to store the MSB info, so significance test is much cheaper. --- include/SPECK2D_INT.h | 1 + include/SPECK2D_INT_ENC.h | 8 +++++++ src/SPECK2D_INT.cpp | 3 +++ src/SPECK2D_INT_ENC.cpp | 45 +++++++++++++++++++++++++++++++-------- 4 files changed, 48 insertions(+), 9 deletions(-) diff --git a/include/SPECK2D_INT.h b/include/SPECK2D_INT.h index a604330f..6ee2946f 100644 --- a/include/SPECK2D_INT.h +++ b/include/SPECK2D_INT.h @@ -44,6 +44,7 @@ class SPECK2D_INT : public SPECK_INT { virtual void m_process_S(size_t idx1, size_t idx2, size_t& counter, bool need_decide) = 0; virtual void m_process_P(size_t idx, size_t& counter, bool need_decide) = 0; virtual void m_process_I(bool need_decide) = 0; + virtual void m_additional_initialization() {}; auto m_partition_S(Set2D) const -> std::array; auto m_partition_I() -> std::array; diff --git a/include/SPECK2D_INT_ENC.h b/include/SPECK2D_INT_ENC.h index 271b0888..ebafb4a4 100644 --- a/include/SPECK2D_INT_ENC.h +++ b/include/SPECK2D_INT_ENC.h @@ -29,9 +29,17 @@ class SPECK2D_INT_ENC final : public SPECK2D_INT { void m_process_S(size_t idx1, size_t idx2, size_t& counter, bool need_decide) final; void m_process_P(size_t idx, size_t& counter, bool need_decide) final; void m_process_I(bool need_decide) final; + void m_additional_initialization() final; + void m_bitplane_init() final; + void m_refinement_extra() final; auto m_decide_S_significance(const Set2D&) const -> bool; auto m_decide_I_significance() const -> bool; + + // `m_msb_buf` stores the MSB bit position of each coefficient, in the same order as + // m_coeff_buf. Significance tests compare entries against `m_msb_threshold`. + std::vector m_msb_buf; + int8_t m_msb_threshold = -1; }; }; // namespace sperr diff --git a/src/SPECK2D_INT.cpp b/src/SPECK2D_INT.cpp index fdb69e3e..5336142b 100644 --- a/src/SPECK2D_INT.cpp +++ b/src/SPECK2D_INT.cpp @@ -212,6 +212,9 @@ void sperr::SPECK2D_INT::m_initialize_lists() m_I.length_x = m_dims[0]; m_I.length_y = m_dims[1]; m_I.part_level = num_of_xforms; + + // Encoder and decoder might have different additional tasks. + m_additional_initialization(); } template class sperr::SPECK2D_INT; diff --git a/src/SPECK2D_INT_ENC.cpp b/src/SPECK2D_INT_ENC.cpp index bc2dee36..b7b82f17 100644 --- a/src/SPECK2D_INT_ENC.cpp +++ b/src/SPECK2D_INT_ENC.cpp @@ -2,6 +2,7 @@ #include #include +#include template void sperr::SPECK2D_INT_ENC::m_process_S(size_t idx1, @@ -31,13 +32,12 @@ void sperr::SPECK2D_INT_ENC::m_process_P(size_t idx, size_t& counter, bool ne bool is_sig = true; if (need_decide) { - is_sig = (m_coeff_buf[idx] >= m_threshold); + is_sig = (m_msb_buf[idx] >= m_msb_threshold); m_bit_buffer.wbit(is_sig); } if (is_sig) { counter++; - m_coeff_buf[idx] -= m_threshold; m_bit_buffer.wbit(m_sign_array.rbit(idx)); m_LSP_new.push_back(idx); m_LIP_mask.wfalse(idx); @@ -65,8 +65,9 @@ auto sperr::SPECK2D_INT_ENC::m_decide_S_significance(const Set2D& set) const assert(!set.is_empty()); for (auto y = set.start_y; y < (set.start_y + set.length_y); y++) { - auto first = m_coeff_buf.data() + y * m_dims[0] + set.start_x; - if (std::any_of(first, first + set.length_x, [th = m_threshold](auto v) { return v >= th; })) + auto first = m_msb_buf.cbegin() + y * m_dims[0] + set.start_x; + if (std::any_of(first, first + set.length_x, + [thld = m_msb_threshold](auto v) { return v >= thld; })) return true; } return false; @@ -79,22 +80,48 @@ auto sperr::SPECK2D_INT_ENC::m_decide_I_significance() const -> bool // It's stored in a contiguous chunk of memory till the buffer end. // assert(m_I.length_x == m_dims[0]); - auto first = m_coeff_buf.data() + size_t{m_I.start_y} * size_t{m_I.length_x}; - auto len = m_coeff_buf.size() - size_t{m_I.start_y} * size_t{m_I.length_x}; - if (std::any_of(first, first + len, [thld = m_threshold](auto v) { return v >= thld; })) + auto first = m_msb_buf.cbegin() + size_t{m_I.start_y} * size_t{m_I.length_x}; + auto len = m_msb_buf.size() - size_t{m_I.start_y} * size_t{m_I.length_x}; + if (std::any_of(first, first + len, + [thld = m_msb_threshold](auto v) { return v >= thld; })) return true; // Second, test the rectangle that's directly to the right of the missing top-left corner. // len = m_dims[0] - m_I.start_x; for (auto y = 0u; y < m_I.start_y; y++) { - first = m_coeff_buf.data() + y * m_dims[0] + m_I.start_x; - if (std::any_of(first, first + len, [thld = m_threshold](auto v) { return v >= thld; })) + first = m_msb_buf.cbegin() + y * m_dims[0] + m_I.start_x; + if (std::any_of(first, first + len, + [thld = m_msb_threshold](auto v) { return v >= thld; })) return true; } return false; } +template +void sperr::SPECK2D_INT_ENC::m_additional_initialization() +{ + const auto len = m_dims[0] * m_dims[1] * m_dims[2]; + m_msb_buf.resize(len); + std::transform(m_coeff_buf.cbegin(), m_coeff_buf.cend(), m_msb_buf.begin(), + [](auto v) { return sperr::msb_position(v); }); +} + +template +void sperr::SPECK2D_INT_ENC::m_bitplane_init() +{ + m_msb_threshold = sperr::msb_position(m_threshold); +} + +template +void sperr::SPECK2D_INT_ENC::m_refinement_extra() +{ + for (auto idx : m_LSP_new) { + assert(m_coeff_buf[idx] >= m_threshold); + m_coeff_buf[idx] -= m_threshold; + } +} + template class sperr::SPECK2D_INT_ENC; template class sperr::SPECK2D_INT_ENC; template class sperr::SPECK2D_INT_ENC; From c9cdb25f83de2adc43bbf85ea3ae3ce394c6048c Mon Sep 17 00:00:00 2001 From: Sam Li Date: Thu, 2 Apr 2026 17:09:53 -0700 Subject: [PATCH 3/3] clang-format, bump version to 0.8.5, and use CLI11 v2.6.2 --- CMakeLists.txt | 4 ++-- src/SPECK2D_INT_ENC.cpp | 6 ++---- 2 files changed, 4 insertions(+), 6 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 186b0587..161fdf4a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.14) -project(SPERR VERSION 0.8.4 LANGUAGES CXX DESCRIPTION "Lossy Scientific Compression with SPERR") +project(SPERR VERSION 0.8.5 LANGUAGES CXX DESCRIPTION "Lossy Scientific Compression with SPERR") if(NOT CMAKE_CXX_STANDARD) set(CMAKE_CXX_STANDARD "20" CACHE STRING "Choose the C++ Standard to use." FORCE) @@ -109,7 +109,7 @@ if( BUILD_CLI_UTILITIES ) set( CLI11_SINGLE_FILE OFF CACHE INTERNAL "Don't use single file CLI11") FetchContent_Declare( cli11 GIT_REPOSITORY https://github.com/CLIUtils/CLI11 - GIT_TAG 6c7b07a878ad834957b98d0f9ce1dbe0cb204fc9 # v2.4.2 + GIT_TAG 37bb6edc5317e99af72ef48405e65d9ca5218861 # v2.6.2 ) FetchContent_MakeAvailable(cli11) diff --git a/src/SPECK2D_INT_ENC.cpp b/src/SPECK2D_INT_ENC.cpp index b7b82f17..1bd18401 100644 --- a/src/SPECK2D_INT_ENC.cpp +++ b/src/SPECK2D_INT_ENC.cpp @@ -82,8 +82,7 @@ auto sperr::SPECK2D_INT_ENC::m_decide_I_significance() const -> bool assert(m_I.length_x == m_dims[0]); auto first = m_msb_buf.cbegin() + size_t{m_I.start_y} * size_t{m_I.length_x}; auto len = m_msb_buf.size() - size_t{m_I.start_y} * size_t{m_I.length_x}; - if (std::any_of(first, first + len, - [thld = m_msb_threshold](auto v) { return v >= thld; })) + if (std::any_of(first, first + len, [thld = m_msb_threshold](auto v) { return v >= thld; })) return true; // Second, test the rectangle that's directly to the right of the missing top-left corner. @@ -91,8 +90,7 @@ auto sperr::SPECK2D_INT_ENC::m_decide_I_significance() const -> bool len = m_dims[0] - m_I.start_x; for (auto y = 0u; y < m_I.start_y; y++) { first = m_msb_buf.cbegin() + y * m_dims[0] + m_I.start_x; - if (std::any_of(first, first + len, - [thld = m_msb_threshold](auto v) { return v >= thld; })) + if (std::any_of(first, first + len, [thld = m_msb_threshold](auto v) { return v >= thld; })) return true; } return false;