Compiler Explorer

Source code

// this is snapshot of project here: https://godbolt.org/z/8f9fMGTKq          
#include <algorithm>
#include <array>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <vector>

#include "dwt53.hpp"
#include "the_matrix.hpp"

#include "utils.hpp"
#include "packmatrix.hpp"
using namespace pack;

the_matrix lenna = {
    {158, 153, 168, 146, 95, 105, 110, 124, 129, 130, 132, 133, 132, 131, 131, 131, 129, 128, 129, 113, 128, 160, 157, 155, 148, 196, 166, 109, 125, 122, 131, 93, },
    {156, 159, 168, 140, 95, 104, 107, 123, 128, 130, 131, 131, 130, 127, 126, 131, 132, 127, 129, 118, 115, 157, 162, 157, 151, 157, 208, 125, 116, 133, 90, 44, },
    {161, 164, 164, 140, 93, 104, 107, 122, 126, 127, 128, 129, 128, 129, 132, 128, 123, 128, 128, 120, 115, 148, 157, 156, 157, 140, 186, 187, 113, 89, 44, 50, },
    {166, 140, 163, 142, 88, 101, 108, 121, 125, 128, 130, 130, 140, 159, 175, 178, 157, 121, 121, 117, 113, 142, 151, 151, 146, 141, 139, 219, 116, 36, 54, 45, },
    {147, 105, 170, 140, 83, 99, 105, 116, 122, 124, 118, 126, 133, 146, 174, 189, 198, 182, 122, 103, 112, 149, 140, 140, 146, 142, 148, 142, 56, 48, 48, 56, },
    {101, 101, 171, 138, 85, 100, 105, 116, 122, 115, 114, 122, 136, 142, 166, 188, 199, 213, 206, 99, 99, 158, 128, 96, 152, 143, 147, 55, 44, 48, 63, 131, },
    {86, 105, 167, 139, 85, 99, 104, 118, 119, 110, 118, 134, 137, 135, 173, 195, 202, 201, 214, 193, 96, 152, 139, 45, 116, 157, 84, 44, 54, 47, 121, 163, },
    {88, 103, 167, 139, 81, 99, 101, 126, 123, 107, 126, 130, 132, 166, 184, 188, 194, 201, 200, 220, 159, 141, 132, 62, 171, 210, 64, 49, 44, 96, 160, 159, },
    {96, 108, 169, 144, 83, 102, 100, 138, 119, 106, 121, 128, 151, 171, 176, 178, 185, 192, 190, 199, 201, 165, 156, 192, 215, 216, 77, 38, 64, 152, 159, 155, },
    {96, 109, 173, 148, 85, 101, 95, 146, 135, 102, 125, 141, 145, 159, 171, 173, 173, 182, 186, 177, 179, 192, 201, 206, 196, 177, 51, 45, 123, 164, 157, 156, },
    {93, 110, 173, 148, 85, 101, 91, 151, 158, 113, 131, 135, 141, 148, 145, 150, 157, 150, 155, 178, 191, 197, 208, 194, 156, 117, 32, 82, 158, 157, 157, 156, },
    {98, 115, 174, 149, 83, 101, 89, 151, 167, 117, 127, 131, 136, 116, 85, 104, 76, 71, 166, 197, 194, 200, 165, 131, 133, 61, 42, 125, 165, 158, 155, 154, },
    {101, 112, 172, 150, 83, 98, 94, 121, 178, 124, 124, 127, 99, 63, 59, 68, 82, 168, 188, 181, 185, 112, 132, 147, 63, 42, 68, 155, 166, 160, 157, 154, },
    {102, 105, 172, 152, 82, 97, 102, 109, 168, 122, 124, 86, 57, 71, 55, 67, 152, 184, 168, 188, 201, 95, 113, 64, 48, 46, 101, 160, 163, 163, 160, 155, },
    {99, 104, 171, 153, 86, 101, 107, 115, 133, 126, 80, 52, 54, 62, 73, 158, 181, 169, 175, 196, 225, 135, 74, 74, 53, 54, 141, 161, 154, 152, 153, 153, },
    {98, 103, 172, 155, 86, 102, 104, 128, 110, 96, 63, 58, 55, 55, 145, 183, 133, 112, 161, 193, 189, 110, 71, 85, 45, 82, 158, 160, 159, 152, 150, 142, },
    {96, 95, 173, 156, 88, 94, 111, 154, 77, 69, 73, 55, 43, 102, 198, 98, 72, 108, 125, 181, 97, 58, 70, 92, 40, 113, 161, 152, 155, 154, 143, 138, },
    {94, 86, 173, 161, 92, 111, 113, 85, 68, 81, 88, 39, 59, 168, 139, 135, 144, 162, 132, 180, 139, 81, 71, 94, 49, 140, 156, 151, 151, 138, 158, 182, },
    {78, 77, 172, 160, 91, 122, 84, 81, 70, 54, 94, 49, 125, 123, 110, 164, 179, 163, 133, 185, 164, 97, 69, 90, 68, 157, 152, 149, 137, 159, 200, 199, },
    {72, 69, 169, 157, 98, 105, 81, 96, 63, 53, 83, 85, 129, 57, 122, 146, 164, 151, 118, 170, 163, 81, 61, 95, 99, 163, 150, 144, 143, 192, 203, 207, },
    {73, 83, 171, 159, 94, 95, 74, 90, 93, 57, 77, 129, 63, 68, 124, 139, 155, 152, 131, 170, 142, 57, 60, 99, 133, 157, 149, 140, 160, 205, 207, 209, },
    {62, 92, 175, 160, 99, 94, 63, 79, 94, 102, 151, 76, 45, 75, 121, 138, 138, 135, 155, 173, 98, 51, 63, 99, 153, 151, 148, 137, 170, 212, 209, 212, },
    {52, 91, 176, 162, 103, 83, 64, 106, 89, 150, 117, 43, 57, 63, 98, 133, 138, 121, 138, 136, 56, 56, 64, 105, 161, 151, 149, 134, 181, 219, 210, 214, },
    {43, 97, 173, 161, 89, 72, 62, 102, 113, 117, 57, 52, 52, 57, 64, 96, 133, 152, 168, 93, 37, 67, 73, 109, 158, 152, 153, 131, 189, 220, 216, 189, },
    {52, 75, 167, 165, 93, 59, 65, 93, 108, 114, 74, 56, 50, 52, 59, 100, 128, 146, 163, 155, 86, 54, 80, 111, 131, 122, 126, 122, 196, 219, 137, 62, },
    {131, 79, 154, 175, 85, 51, 73, 75, 84, 119, 102, 58, 52, 57, 69, 115, 137, 140, 161, 196, 210, 136, 55, 117, 149, 131, 112, 113, 210, 158, 69, 92, },
    {121, 113, 151, 177, 73, 45, 67, 55, 72, 110, 115, 70, 50, 60, 85, 113, 140, 143, 158, 172, 194, 228, 91, 110, 153, 140, 129, 172, 215, 107, 99, 102, },
    {98, 116, 149, 175, 66, 46, 70, 50, 57, 106, 114, 78, 54, 75, 104, 118, 141, 139, 151, 167, 188, 221, 161, 121, 154, 131, 139, 221, 176, 98, 102, 96, },
    {82, 123, 147, 176, 72, 48, 59, 52, 64, 87, 88, 60, 48, 102, 104, 132, 142, 137, 146, 160, 182, 207, 204, 122, 128, 141, 137, 200, 145, 86, 97, 91, },
    {91, 137, 146, 175, 72, 45, 59, 50, 57, 98, 73, 38, 67, 99, 103, 143, 142, 138, 140, 153, 172, 197, 219, 122, 81, 68, 154, 186, 82, 86, 96, 79, },
    {108, 146, 151, 171, 61, 49, 56, 68, 56, 89, 64, 39, 80, 98, 130, 138, 142, 144, 140, 146, 162, 184, 217, 159, 94, 86, 144, 101, 74, 85, 100, 69, },
    {73, 164, 159, 162, 51, 58, 65, 79, 75, 72, 63, 69, 113, 128, 131, 135, 140, 145, 146, 147, 158, 175, 206, 182, 92, 123, 114, 96, 83, 100, 89, 59, },
};

int main() {
    int max_levels = 5;
    dwt2d::Wavelet wavelet = dwt2d::dwt53;
    //auto data = make_envelope(32,32,1);
    // cubicBlur3x3(data);
    // cubicBlur3x3(data);
    //auto data = make_gradient(32,32,0,11,11,22);
    auto data = lenna;
    //auto data = make_sky(32,32);
    //auto data = make_random(32);
    // cubicBlur3x3(data);
    // cubicBlur3x3(data);

std::cout << "original: pw=" << matrix_energy(data) << raster(data);
    auto data_comp = compress(data);
    std::cout << "packed by huffman size = " << data_comp.size() << std::endl;
    // auto data_decomp = huffman::decompress(data_comp);
    // std::cout <<  data_decomp;

// auto test = make_matrix(10,10, [](int x,int y,int&v){ v = x + y*10; });
    //  auto all_img = make_block( data );
    //   quantization_block(data,all_img, true);
    //   quantization_block(data,all_img, false);
    //  std::cout << "quantized:" << raster(data);
    //  std::cout << "psnr=" << psnr(lenna, data) << std::endl;

dwt2d::Transform codec;
    codec.prepare_transform(max_levels, wavelet, data);
    auto& haar_data = codec.forward();

// std::cout << raster(haar_data);
    std::cout << "transformed: pw=" << matrix_energy(haar_data) << haar_data;

auto haar_data_comp = compress(haar_data);
    decompress(haar_data_comp);
    std::cout << "packed by huffman size = " << haar_data_comp.size() << std::endl;
    
    codec.vq_forward();
    std::cout << "quantized: pw=" << matrix_energy(haar_data) << haar_data;

auto haar_data_vq_comp = compress(haar_data);
    decompress(haar_data_vq_comp);
    std::cout << "packed by huffman size = " << haar_data_vq_comp.size() << std::endl;

codec.vq_inverse();
    std::cout << "quantized-invert: pw=" << matrix_energy(haar_data) << haar_data;

// std::cout << haar_data;
    auto& reconstructed = codec.inverse();

std::cout << "psnr=" << psnr(data, reconstructed) << raster(reconstructed);

auto depacked = decompress(data_comp);
    std::cout << "data unpacked:" << raster(depacked);
}

Source code

#pragma once
#include <vector>
#include <iostream>
#include <cstdint>
#include <algorithm>
#include <vector>
#include "the_matrix.hpp"
#include "huffman.hpp"
#include "utils.hpp"

namespace pack {

inline std::vector<uint8_t> compress(const the_matrix& matrix) {
    auto data = flatten(matrix);
    Huffman<int> henc;
    BitWriter enc_w;
    henc.buildHuffmanTree(data.begin(), data.end());
    henc.encodeHuffmanTree(enc_w);
        
    std::cerr << "haffmantable size = " << enc_w.size() << std::endl;

enc_w.writeBits(16, matrix[0].size());
    enc_w.writeBits(16, matrix.size());    
    henc.encode(enc_w,data.begin(),data.end());
    enc_w.flush();

return std::vector<uint8_t>(enc_w.data(), enc_w.data() + enc_w.size());
}

inline the_matrix decompress(const std::vector<uint8_t>& compressed) {
    Huffman<int> hdec;
    BitReader enc_r(compressed.data(),compressed.size());
    hdec.decodeHuffmanTree(enc_r);   
    return make_matrix( enc_r.readBits(16), enc_r.readBits(16),
    [&](int &v) {
        v = hdec.decode(enc_r);
    }
    );
}

}  // namespace pack

Source code

//https://raw.githubusercontent.com/vovach777/BitStream/main/bitstream.hpp
#pragma once
#include <cstdint>
#include <iostream>
#include <ostream>
#include <vector>
#include <stdexcept>

class BitWriter {

inline void put_bits(int n, uint32_t value) {
    if (n == 0)
        throw std::out_of_range("n == 0, can not write 0 bits");
    if (   ~static_cast<uint32_t>((1<<n)-1) & value) {
        std::cerr << std::endl << "bits = " << n << ", value = " << value << std::endl;
        throw std::out_of_range("value out of n-bits !!!");
    }

auto bit_left = bit_left_;
    auto bit_buf  = bit_buf_;
    //allways reserve
    if (vec.size() <= size_) {
        vec.resize(size_+1);
    }

if (n < bit_left) {
        bit_buf     = (bit_buf << n) | value;
        bit_left   -= n;
    } else {
        bit_buf   <<= bit_left;
        bit_buf    |= value >> (n - bit_left);
        vec[size_++] = __builtin_bswap32(bit_buf);
        bit_left   += 32 - n;
        bit_buf     = value;
    }
    bit_buf_  = bit_buf;
    bit_left_ = bit_left;
    }

public:
    inline void writeBits(int n, uint32_t value ) { put_bits(n, value); }
    inline void writeBit(bool value) { put_bits(1,value);}
    inline void writeBit0() { writeBit(true);}
    inline void writeBit1() { writeBit(false);}
    auto data() const { 
        return  reinterpret_cast<const uint8_t*>( vec.data()); 
    }
    auto data() { 
        return  reinterpret_cast<uint8_t*>( vec.data()); 
    }
    auto size_in_bits() const { return size_*32 + (32-bit_left_); }
    auto size() const {  
       return size_in_bits() + 7 >> 3;
        }
  
    void flush() {
        if (bit_left_ == 32)
            return;
   
        vec[size_] = 0;
        auto byte_by_byte = reinterpret_cast<uint8_t*>( &vec[size_] );     
        auto bit_left = bit_left_;
        auto bit_buf = bit_buf_;
        if (bit_left < 32)
            bit_buf <<= bit_left;
        while (bit_left < 32) {
            *byte_by_byte++ = bit_buf >> (32 - 8);
            bit_buf  <<= 8;
            bit_left  += 8;
        }
    }

// Метод для вывода данных (для проверки)
    void printBits( std::ostream& o) const {
        auto count = size_in_bits();
        for (auto it = data(), end=data()+size(); it != end; ++it ) {
            auto byte = *it;
            for (int i = 7; i >= 0; --i) {
                  if (count-- ==  0)
                    return;
                o << ((byte >> i) & 1);
            }
            o << " ";
        }
        //o << "\n";
    }

private:
    std::vector<uint32_t> vec;
    uint32_t bit_buf_ = 0;  // буфер для хранения битов
    int bit_left_ = 32;  // счетчик битов в буфере
    std::vector<uint32_t>::size_type size_=0; //current vector tail
};

class BitReader {

inline unsigned int get_bits1() {
        size_t _index = index;
        if (index >= size_in_bits) {
            throw std::out_of_range("(index >= size_in_bits)!!!");
        }
        uint8_t result = buffer[_index >> 3];
        result <<= _index & 7;
        result >>= 8 - 1;
        _index++;
        index = _index;
        return result;
    }

/**
     * Read 1-25 bits.
     */
    inline unsigned int get_bits(int n) {
        // assert(n>0 && n<=25);
        if (n == 0) return 0;
        if (index+n > size_in_bits) {
            throw std::out_of_range("(index+n > size_in_bits)!!!");
        }       
        union unaligned_32 {
            uint32_t l;
        } __attribute__((packed)) __attribute__((may_alias));
        unsigned int tmp;
        unsigned int re_index = index;
        unsigned int re_cache;
        re_cache =
            __builtin_bswap32(
                ((const union unaligned_32 *)(buffer + (re_index >> 3)))->l)
            << (re_index & 7);
        tmp = (((uint32_t)(re_cache)) >> (32 - n));
        re_index += n;
        index = re_index;
        return tmp;
    }

/**
     * Read 0-32 bits.
     */

inline unsigned int get_bits_long(int n) {
        // assert(n>=0 && n<=32);
        if (!n) {
            return 0;
        if (n==1)
            return get_bits1();

} else if (n <= 25) {
            return get_bits(n);
        } else {
            unsigned ret = get_bits(16) << (n - 16);
            return ret | get_bits(n - 16);
        }
    }

public:
    BitReader(const uint8_t *buffer_, size_t buffer_size, size_t limit_bits=std::numeric_limits<size_t>::max()) : size_in_bits(std::min(buffer_size << 3,limit_bits)) {
        buffer_begin = buffer = buffer_;
        buffer_end = buffer_ + buffer_size;
        index = 0;
    }
    BitReader() = default;
    inline bool readBit() {
        return get_bits1();
    }
    inline auto readBits(int n) {
        return get_bits_long(n);
    }
    // BitReader& operator=(BitReader&& other) {

// }
    // Метод для вывода данных (для проверки)
    void printBits(std::ostream & o) const {
        auto count  = size_in_bits;
        for (auto it = buffer_begin; it != buffer_end; ++it) {
            auto byte = *it;
            for (int i = 7; i >= 0; --i) {
                if (count-- == 0)
                    return;
                o << ((byte >> i) & 1);
            }
            o << " ";
        }
       // o << "\n";
    }

size_t get_index() { return index; }
    //size_t size_in_bits() { return size_in_bits; }
    size_t bit_left() { return size_in_bits - index; }

private:
    const uint8_t *buffer_begin;
    const uint8_t *buffer;
    const uint8_t *buffer_end;
    size_t index;
    size_t size_in_bits;
};

inline std::ostream& operator<<(std::ostream& o, BitWriter const & a) {
    a.printBits(o);
    return o;
}

inline std::ostream& operator<<(std::ostream& o, BitReader const& a) {
    a.printBits(o);
    return o;
}

Source code

#pragma once
/*
  * variants *
        ** L - variants **
        A: xLH
        B: HLx
        D: HLH

** H - variants **
        C: LHL
        E: LHx

** ODD N=7 **
        LHLHLHL => H-variants:  LHL (C);   L-variants:  xLH (A), HLH (D), HLx
  (B)

** EVEN N=8 **
        LHLHLHLH => H-variants: LHL (C), LHx (E); L-variants: xLH (A), HLH (D)

** Formulas (forward) **
        H, L
        C: H = H - (L0 + L1) / 2
        E: H = H - L0
        D: L = L + (H0 + H1 + 2) / 4
        B: L = L + (H0 + 1) / 2
        A: L = L + (H1 + 1) / 2
        ** Formulas (inverse) **
        L, H
        D: L = L - (H0 + H1 + 2) / 4
        B: L = L - (H0 + 1) / 2
        A: L = L - (H1 + 1) / 2
        C: H = H + (L0 + L1) / 2
        E: H = H + L0
*/

#ifdef __cplusplus
#include "the_matrix.hpp"
#endif

static int sizeof_L(int x) { return (x + 1) >> 1; }

static int sizeof_H(int x) { return x >> 1; }

static void dwt53(const int *x, int size, int *L, int *H) {
    if (size & 1) {
        // C-variant
        int *H_p = H;
        const int *x_p = x + 1;
        int s = size / 2;
        while (s--) {
            *H_p++ = x_p[0] - (x_p[-1] + x_p[1]) / 2;
            x_p += 2;
        }
        int *L_p = L;
        // A-variant
        *L_p++ = x[0] + (H[0] + 1) / 2;
        // D-variant
        s = (size + 1) / 2 - 2;
        x_p = x + 2;
        H_p = H;
        if (s == 0)
            H_p++;
        else
            while (s--) {
                *L_p++ = x_p[0] + (H_p[0] + H_p[1] + 2) / 4;
                H_p += 1;
                x_p += 2;
            }
        // B-variant
        *L_p++ = x_p[0] + (H_p[-1] + 1) / 2;

} else {
        // C-variant
        int s = size / 2 - 1;
        int *H_p = H;
        const int *x_p = x + 1;
        while (s--) {
            *H_p++ = x_p[0] - (x_p[-1] + x_p[1]) / 2;
            x_p += 2;
        }
        // E-variant
        *H_p = x_p[0] - x_p[-1];
        s = size / 2 - 1;
        int *L_p = L;
        H_p = H;
        // A-variant
        *L_p++ = x[0] + (H_p[0] + 1) / 2;
        // D-variant
        x_p = x + 2;
        while (s--) {
            *L_p++ = x_p[0] + (H_p[0] + H_p[1] + 2) / 4;
            H_p += 1;
            x_p += 2;
        }
    }
}

static void idwt53(int *x, int size, const int *L, const int *H) {
    if (size & 1) {
        const int *L_p = L;
        // A-variant
        //*L_p++ = x[0]+(H[0] + 1)/2;
        x[0] = *L_p++ - (H[0] + 1) / 2;
        // D-variant
        int s = (size + 1) / 2 - 2;
        int *x_p = x + 2;
        const int *H_p = H;
        if (s == 0)
            H_p++;
        else
            while (s--) {
                //*L_p++ = x_p[0] + (H_p[0] + H_p[1] + 2) / 4;
                x_p[0] = *L_p++ - (H_p[0] + H_p[1] + 2) / 4;
                H_p += 1;
                x_p += 2;
            }
        // B-variant
        //*L_p++ = x_p[0]+(H_p[-1] + 1) /2;
        x_p[0] = *L_p++ - (H_p[-1] + 1) / 2;

// C-variant
        H_p = H;
        x_p = x + 1;
        s = size / 2;
        while (s--) {
            //*H_p++ = x_p[0] - (x_p[-1] + x_p[1])/2;
            x_p[0] = *H_p++ + (x_p[-1] + x_p[1]) / 2;
            x_p += 2;
        }
    } else {
        // LLLLLL
        //////
        int s = size / 2 - 1;
        const int *L_p = L;
        const int *H_p = H;
        // A-variant
        x[0] = *L_p++ - (H_p[0] + 1) / 2;
        // D-variant
        int *x_p = x + 2;
        while (s--) {
            x_p[0] = *L_p++ - (H_p[0] + H_p[1] + 2) / 4;
            H_p += 1;
            x_p += 2;
        }
        // HHHHHHH
        // C-variant
        s = size / 2 - 1;
        H_p = H;
        x_p = x + 1;
        while (s--) {
            x_p[0] = *H_p++ + (x_p[-1] + x_p[1]) / 2;
            x_p += 2;
        }
        // E-variant
        x_p[0] = *H_p + x_p[-1];
    }
}

static void haar(const int *x, int size, int *L, int *H) {
    int *H_p = H;
    int *L_p = L;
    const int *x_p = x;
    int s = size / 2;
    while (s--) {
        int a = x_p[0];
        int b = x_p[1];
        b -= a;
        a += b / 2;
        *L_p++ = a;
        *H_p++ = b;
        x_p += 2;
    }
    if (size & 1) {
        *L_p++ = x_p[0];
    }
}

static void ihaar(int *x, int size, const int *L, const int *H) {
    const int *H_p = H;
    const int *L_p = L;
    int *x_p = x;
    int s = size / 2;
    while (s--) {
        int a = *L_p++;
        int b = *H_p++;
        a -= b / 2;
        b += a;
        *x_p++ = a;
        *x_p++ = b;
    }
    if (size & 1) {
        *x_p = *L_p;
    }
}

#ifdef THE_MATRIX_DEFINED
#include "quantization.hpp"

namespace dwt2d {

using DWTFunction = decltype(dwt53) *;
using IDWTFunction = decltype(idwt53) *;

enum Wavelet { unchanged, haar, dwt53 };

enum BlockType { HL, LH, HH };

struct Band {
    int nb;
    int width;
    int height;
    int widthLL;
    int heightLL;
    Block blocks[3];
    Band(int num, int &width, int &height)
        : nb(num), width(width), height(height) {
        const auto L_width = sizeof_L(width);
        const auto H_width = sizeof_H(width);
        const auto L_height = sizeof_L(height);
        const auto H_height = sizeof_H(height);
        widthLL = L_width;
        heightLL = L_height;

blocks[HL].x.offs = L_width;
        blocks[HL].x.size = H_width;
        blocks[HL].y.offs = 0;
        blocks[HL].y.size = L_height;

blocks[LH].x.offs = 0;
        blocks[LH].x.size = L_width;
        blocks[LH].y.offs = L_height;
        blocks[LH].y.size = H_height;

blocks[HH].x.offs = L_width;
        blocks[HH].x.size = H_width;
        blocks[HH].y.offs = L_height;
        blocks[HH].y.size = H_height;

width = L_width;
        height = L_height;
    }
};

static std::vector<dwt2d::Band> calc_geometry(int width, int height,
                                              int levels) {
    std::vector<Band> result;
    int num = 1;
    while (levels-- && width > 1 && height > 1) {
        auto &lb = result.emplace_back(num++, width, height);
    }
    return result;
}

template <typename F>
static void process_levels(the_matrix &a, int levels, F f) {
    int width = a[0].size();
    int height = a.size();

for (int level = 1; level <= levels && width > 1 && height > 1; ++level) {
        dwt2d::Band band(level, width, height);
        f(band);
    }
}

class Transform {
    the_matrix data_;
    DWTFunction forward_ = &::dwt53;
    IDWTFunction inverse_ = &::idwt53;
    int levels_ = 1;
    // std::vector<int> x_sizes_;
    // std::vector<int> y_sizes_;
    std::vector<int> LH_storage_;
    std::vector<dwt2d::Band> geometry_;
    std::vector<int> hist_;
    Codebook codebook_;

public:
   Transform() {

}
    void set_transform(DWTFunction forward, IDWTFunction inverse) {
        forward_ = forward;
        inverse_ = inverse;
    }
    void set_transform(Wavelet transform) {
        if (transform == haar) {
            forward_ = &::haar;
            inverse_ = &::ihaar;
        } else if (transform == dwt53) {
            forward_ = &::dwt53;
            inverse_ = &::idwt53;
        }
    }

the_matrix &get_data() { return data_; }

// static std::vector<int> calc_L_sizes( int size, int levels)
    // {
    //     std::vector<int> sizes(levels);
    //     sizes.clear();
    //     for (int level = 1; level <= levels && (size > 1); level++, size =
    //     sizeof_L(size))
    //     {
    //         sizes.push_back( size );
    //     }
    //     return sizes;
    // }

void prepare_transform(int levels, Wavelet transform,
                           the_matrix const &data) {
        data_ = data;
        set_transform(transform);
        // x_sizes_ = calc_L_sizes(data_[0].size());
        // y_sizes_ = calc_L_sizes(data_.size());
        LH_storage_.resize(std::max(data[0].size(), data.size()));
        geometry_ = calc_geometry(data[0].size(), data.size(), levels);
        levels_ = geometry_.size();
        hist_.clear();
        codebook_.clear();
    }

the_matrix &forward() {
        for (auto &row : data_) {
            for (int i = 0; i < geometry_.size(); i++) {
                const auto size = geometry_[i].width;
                forward_(row.data(), size, LH_storage_.data(),
                         LH_storage_.data() + sizeof_L(size));
                std::copy(LH_storage_.begin(), LH_storage_.begin() + size,
                          row.begin());
            }
        }
        for (int x = 0; x < data_[0].size(); x++) {
            auto col = get_col(data_, x);
            for (int i = 0; i < geometry_.size(); i++) {
                const auto size = geometry_[i].height;
                forward_(col.data(), size, LH_storage_.data(),
                         LH_storage_.data() + sizeof_L(size));
                std::copy(LH_storage_.begin(), LH_storage_.begin() + size,
                          col.begin());
            }
            set_col(data_, x, col);
        }
        return data_;
    }

the_matrix &inverse() {
        for (int x = 0; x < data_[0].size(); x++) {
            auto col = get_col(data_, x);
            for (int i = geometry_.size() - 1; i >= 0; --i) {
                const auto size = geometry_[i].height;
                std::copy(col.begin(), col.begin() + size, LH_storage_.begin());
                inverse_(col.data(), size, LH_storage_.data(),
                         LH_storage_.data() + sizeof_L(size));
            }
            set_col(data_, x, col);
        }

for (auto &row : data_) {
            for (int i = geometry_.size() - 1; i >= 0; --i) {
                const auto size = geometry_[i].width;
                std::copy(row.begin(), row.begin() + size, LH_storage_.begin());
                inverse_(row.data(), size, LH_storage_.data(),
                         LH_storage_.data() + sizeof_L(size));
            }
        }
        return data_;
    }

void vq_forward() {
        hist_.clear();
        hist_.reserve(1024);
        // histogramm
        for (int level = 0; level < levels_; ++level) {
            for (int i = 0; i < 3; i++)
                process_matrix(data_, geometry_[level].blocks[i],
                               [this](int v) {
                                   int value = v;
                                   if (value < 0) value = -value;
                                   if (value >= hist_.size()) {
                                       hist_.resize(value + 1);
                                   }
                                   ++hist_[value];
                               });
        };

codebook_ = Codebook(hist_);
        std::cerr << "codebook=" << codebook_.get_codebook();

for (int level = 0; level < levels_; ++level) {
            for (int i = 0; i < 3; i++)
                process_matrix(
                    data_, geometry_[level].blocks[i],
                    [this](int &v) { v = codebook_.get_codeword(v);});
        }
    }

void vq_inverse() {
        if (codebook_.size()) {
            for (int level = 0; level < levels_; ++level) {
                for (int i = 0; i < 3; i++)
                    process_matrix(
                        data_, geometry_[level].blocks[i],
                        [this](int &v) { v = codebook_.get_value(v); });
            }
        }
    }

void vq_inverse(Codebook codebook) {
        codebook_ = std::move(codebook);
        vq_inverse();
    }
};

}  // namespace dwt2d

#endif

Source code

//https://raw.githubusercontent.com/vovach777/huffman_advanced/main/utils.hpp
#pragma once

#include <algorithm>

template <typename T>
inline T square(const T v) {
   return v*v;
}

template <typename T>
inline T median(T a, T b, T c) {
   if (a<b) {
      if (b<c)
         return b;
      return  (a < c) ? c : a;
   } else {
      if (b>=c)
         return b;
      return a < c ? a : c;
   }
}

inline int r_shift( int num, int shift) {
    
    return (num < 0 ? -(-num >> shift) : num >> shift);
    
}

template <typename Iterator>
Iterator find_nerest(Iterator begin, Iterator end, typename std::iterator_traits<Iterator>::value_type val)
{
    if (begin == end)
        return end;
    auto it = std::lower_bound(begin, end, val);

if (it == end ) {
        return std::prev(end);
    };
    if (it == begin)
        return it;
    if (*it == val)
        return it;

const auto lower_distance = std::abs(*std::prev(it) - val);
    const auto upper_distance = std::abs(*it - val);
    return lower_distance < upper_distance ? std::prev(it) : it;
}

// Function to convert a coefficient into a "category/index" pair in JPEG
    // format
   inline int symbol_to_catindex(int coeff) {
        // If the coefficient is negative, we change its sign to positive
        uint32_t positive = (coeff < 0) ? -coeff : coeff;
        if (positive == 0) return 0;

// Calculate the category of the coefficient as the number of bits
        // needed to represent its absolute value
        int cat = 32 - __builtin_clz(positive);

// Calculate the minimum and maximum value in the given category
        int minValue = (1 << (cat - 1));
        int maxValue = (1 << cat) - 1;

// Calculate the index of the coefficient within the given category
        auto index = (coeff < 0) ? maxValue + coeff
                                 : coeff - minValue + (maxValue - minValue + 1);

// Return the "category/index" pair in JPEG format
        return cat | index << 4;
    }

// Function to recover the coefficient from a "category/index" pair in JPEG
    // format
    inline int catindex_to_symbol(int pair) {
        if (pair == 0) return 0;

// Extract the category and index from the "category/index" pair
        int cat = pair & 0xf;
        int index = pair >> 4;

// Calculate the minimum and maximum value in the given category
        int minValue = (1 << (cat - 1));
        int maxValue = (1 << cat) - 1;

// Calculate the index of the positive value within the given category
        auto positive_index = 1 << (cat - 1);

// Recover the original coefficient
        return (index < positive_index) ? -maxValue + index
                                        : index - positive_index + minValue;
    }

// Function to recover the coefficient from a "category/index" pair in JPEG
    // format
    inline int catindex_to_symbol(int cat, int index) {
        if (cat == 0) return 0;

// Calculate the minimum and maximum value in the given category
        int minValue = (1 << (cat - 1));
        int maxValue = (1 << cat) - 1;

// Calculate the index of the positive value within the given category
        auto positive_index = 1 << (cat - 1);

// Recover the original coefficient
        return (index < positive_index) ? -maxValue + index
                                        : index - positive_index + minValue;
    }

Source code

#ifndef THE_MATRIX_DEFINED
   #define THE_MATRIX_DEFINED

#define _USE_MATH_DEFINES
#include <cmath>
#include <vector>
#include <algorithm>
#include <ostream>
#include <iomanip>
#include <cassert>
#include <cstdint>

struct Index {
    int offs;
    int size;
};

struct Block {
    Index x;
    Index y;
};

using the_matrix = std::vector<std::vector<int>>;
static bool raster_flag = false;

static int make_positive(the_matrix & a);
static inline int getValue(const the_matrix& matrix, int row, int col);

template <typename MatrixType>
decltype(auto) raster(MatrixType && img) {
    raster_flag = true;
    return std::forward<MatrixType>(img);
}

struct Range {
    int min;
    int max;
};

inline Block make_block(the_matrix const & matrix){
    //return Block{{0,matrix[0].size()}, {0,matrix.size()}};
    Block result;
    result.x.offs = 0;
    result.x.size = matrix[0].size();
    result.y.offs = 0;
    result.y.size = matrix.size();
    return result;
}

inline Range get_range(the_matrix const &matrix)
{
    int minVal = std::numeric_limits<int>::max();
    int maxVal = std::numeric_limits<int>::min();

for (const auto& row : matrix) {
        for (const auto val : row) {
            minVal = std::min(minVal, val);
            maxVal = std::max(maxVal, val);
        }
    }
    return Range{minVal,maxVal};
}

inline Range get_range(the_matrix const &matrix, Block const& block)
{
    int minVal = std::numeric_limits<int>::max();
    int maxVal = std::numeric_limits<int>::min();
    
    for (int j=0; j < block.y.size; j++)
    for (int i=0; i < block.x.size; i++)
    {
        auto val = getValue(matrix, block.x.offs + i, block.y.offs + j);
        minVal = std::min(minVal, val);
        maxVal = std::max(maxVal, val);        
    }
    return Range{minVal,maxVal};
}

inline void raster(std::ostream& o,the_matrix const& img)
{
    static char pixels[] = {' ', ':', ';' ,'x','X','@'};
    auto range = get_range(img);

for (const auto & row : img) {
      for (const auto v : row)
      {
         auto pixel = (range.max - range.min == 0) ? pixels[0] : pixels[ std::clamp<int>(  (v - range.min) * sizeof(pixels) / (range.max - range.min),0,sizeof(pixels)-1) ];
         o <<  pixel << pixel << pixel ;
      }
        o <<std::endl;
      }
}

inline std::ostream& raster(std::ostream& os) {
    raster_flag = true;
    return os;
}

inline std::ostream& operator << (std::ostream& o,  std::vector<int> const & a)
{   
    for (auto v : a)
        o << std::setw(3) << v << " ";
    o << std::endl;
    return o;
}

inline std::ostream& operator << (std::ostream& o,  the_matrix const & a)
{
     o << std::endl;
    if (raster_flag) {
        raster_flag = false;
        raster(o, a);
    }
    else {
        for (auto const & row : a)
        {
            o << row;
        }
    }
    o << std::endl;
    return o;
}

inline void transpose(the_matrix& matrix) {
    const auto rows = matrix.size();
    const auto cols = matrix[0].size();

if (rows != cols) {
        // Матрица не квадратная, выполняем перемещение
        the_matrix transposed(cols, std::vector<int>(rows));

for (int i = 0; i < rows; ++i) {
            for (int j = 0; j < cols; ++j) {
                transposed[j][i] = matrix[i][j];
            }
        }

std::swap(matrix,transposed);
    } else {
        // Матрица квадратная, выполняем обмен элементами
        for (int i = 0; i < rows; ++i) {
            for (int j = 0; j < i; ++j) {
                std::swap(matrix[i][j], matrix[j][i]);
            }
        }
    }
}

inline std::vector<int> get_col(the_matrix const &matrix, int x) {
    std::vector<int> col(matrix.size());
    for (int i=0; i< col.size(); i++ )
        col[i] = matrix[i][x];
    return col;
}

inline void set_col(the_matrix  &matrix, int x,  std::vector<int> const & col) {
    assert( col.size() == matrix.size() );
    
    for (int i=0; i< col.size(); i++ )
        matrix[i][x] = col[i];
    
}

inline void setValue( the_matrix& matrix, int row, int col, int value) {
    int height = matrix.size();
    int width = matrix[0].size();

if (row >= 0 && row < height && col >= 0 && col < width)
        matrix[row][col] = value;
}

inline void drawLine(the_matrix& matrix, int x1, int y1, int x2, int y2, int value = 1) {
    int width = matrix[0].size();
    int height = matrix.size();

int dx = x2-x1;
    int dy = y2-y1;
    double m = double(dy)/dx;
    double xin,yin;
    int tot;
    if(dx==0)
    {
        if(y1>y2) std::swap(y1,y2);
        for(int i=y1;i<=y2;i++) setValue(matrix,i,x1,value);
        return ;
    }
    else if(dy==0)
    {
        if(x1>x2) std::swap(x1,x2);
        for(int i=x1;i<=x2;i++) setValue(matrix,y1,i,value);
        return ;
    }

if(abs(dx)>=abs(dy))
    {
        if(x1>x2)
        {
            std::swap(x1,x2);
            std::swap(y1,y2);
        }
        xin = 1;
        yin = m;
        tot = x2-x1+1;
    }
    else
    {
        if(y1>y2)
        {
            std::swap(x1,x2);
            std::swap(y1,y2);
        }
        yin = 1;
        xin = 1.0/m;
        tot = y2-y1 + 1;

}
    double x = x1;
    double y = y1;
    for(int i=1;i<=tot;i++)
    {
        setValue(matrix,round(y),round(x),value);
        x += xin;
        y += yin;
    }

}

inline int getValue(const the_matrix& matrix, int row, int col) {
    int height = matrix.size();
    int width = matrix[0].size();

// Удерживаем координаты в границах матрицы
    row = std::max(0, std::min(height - 1, row));
    col = std::max(0, std::min(width - 1, col));

return matrix[row][col];
}

inline void cubicBlur3x3(the_matrix& matrix) {
    int height = matrix.size();
    int width = matrix[0].size();

// Копирование исходной матрицы для сохранения оригинальных значений
    the_matrix originalMatrix = matrix;

// Применение кубического размытия
    for (int y = 0; y < height; y++) {
        for (int x = 0; x < width; x++) {
            int sum = 0;
            for (int yy = y-1; yy <= y+1; yy++) {
                for (int xx = x-1; xx <= x+1; xx++) {
                    sum += getValue(originalMatrix, yy, xx);
                }
            }
            matrix[y][x] = (sum+5) / 9;
        }
    }
}

inline the_matrix make_radial(int size) {
   the_matrix data(size,std::vector<int>(size));
  for (int y=0; y<size; y++)
  for (int x=0; x<size; x++)
    {
        double xx =  x * 2.1 - size;
        double yy =  y * 2.1 - size;
        double r = sqrt( xx*xx+yy*yy );
        //data[i] = r > 8.0 ? 0 : (int)r * 16 / 8;
        data[y][x] = (int)r*256/size/2.1;
        //data[i] = i % 8;
    }
    return data;
}

inline the_matrix make_random(int size) {
    the_matrix data(size,std::vector<int>(size));
    for (int y=0; y<size; y++)
    for (int x=0; x<size; x++)
    {

data[y][x] = rand() % 256;
        //data[i] = i % 8;
    }
    return data;
}

inline the_matrix make_monotonic(int size, int value)
{
    return the_matrix(size, std::vector<int>(size,value));
}

inline the_matrix make_gradient(int size) {
   the_matrix data(size,std::vector<int>(size));
  for (int y=0; y<size; y++)
  for (int x=0; x<size; x++)
    {
        data[y][x] = 100 + (x+y)*100/(size*2);
        //data[i] = i % 8;
    }
    return data;
}

inline std::vector<double> make_sin(int N, int freq = 1, double phase_in_T=0) {

std::vector<double> vec;

int phase_in_N = (int)round( N * fmod( phase_in_T, 1) / freq) % N;
    for (int i=0; i<N; i++)
        vec.push_back( std::sin( 2 * M_PI * freq / N * ( (i+phase_in_N) % N)  ) );
    return vec;
}

inline double interpolate(
    const double bottomLeft,
    const double topLeft,
    const double bottomRight,
    const double topRight,
    const double dx,
    const double dy
) {
    return
        bottomLeft * (1 - dx) * (1 - dy) +
        topLeft    * (1 - dx) * dy +
        bottomRight * dx * (1 - dy) +
        topRight   * dx * dy;
}

inline the_matrix make_gradient(int width, int height, int topLeft, int topRight, int bottomLeft, int bottomRight) {
    the_matrix a( height, std::vector<int>(width,0));

for (int y = 0; y < height; ++y) {
        for (int x = 0; x < width; ++x) {
            double dx = static_cast<double>(x) / (width - 1);
            double dy = static_cast<double>(y) / (height - 1);
            a[y][x] = interpolate(topLeft, bottomLeft, topRight,bottomRight, dx, dy);
        }
    }
    return a;
}

inline int make_positive(the_matrix & a)
{
    int min = 0;
    for (const auto & row: a)
    for (const auto v : row)
       min = std::min(v, min);
    if (min < 0) {
    for (auto & row: a)
    for (auto & v : row)
       v -= min;
    }
    return min;
}

inline the_matrix make_sin2d(int width, int height, double freq, int time) {
    the_matrix gradient = the_matrix(height, std::vector<int>(width,0));
    double freq_x = freq + sin(time/16.0);
    double freq_y = freq_x / 2;
    double amp_y = 12;
    for (int y=0, time_y=time; y < height; ++y, ++time_y, amp_y *= 0.98 ) {
        double amp_x = 64;

int y_ = static_cast<int>( std::sin( 2 * M_PI * freq_y / height * (time_y % height)) * amp_y );
        for (int t = time+width, x=0; x < width; --t,++x, amp_x *= 0.98  )
        {
            ////gradient[y][x]  = y_ - static_cast<int>( y * 0.2 );
            gradient[y][x]  = y_  +  static_cast<int>( std::sin( 2 * M_PI * freq_x / width * (t % width)) * amp_x ) + 128;
        }
    }
   // make_positive(gradient);
    return gradient;
}

inline the_matrix make_envelope(int width, int height, int value)
{
    the_matrix result(height, std::vector<int>(width, 0));
    drawLine(result, 0, 0, width-1, 0, value);
    drawLine(result, width-1, 0, width-1, height-1, value);
    drawLine(result, width-1, height-1, 0, height-1, value);
    drawLine(result, 0, height-1, 0, 0, value);
    drawLine(result, 0, 0, width-1, height-1, value);
    drawLine(result, 0, height-1, width-1, 0, value);
    return result;
}

inline int psnr(const the_matrix & originalImage, const the_matrix& compressedImage)
{
    // Расчет суммы квадратов разностей пикселей
    auto mse = 0ULL;
    auto originalRange = get_range(originalImage);
    int count = originalImage.size() * originalImage[0].size();
    if (count == 0)
        return 0;
    for (int y=0,max_y=std::min(originalImage.size(), compressedImage.size()); y<max_y; y++)
    for (int x=0,max_x=std::min(originalImage[y].size(), compressedImage[y].size()); x<max_x; x++)
    {
        auto diff = originalImage[y][x]+originalRange.min - (compressedImage[y][x] + originalRange.min);
        mse += diff * diff;
    }

// Расчет среднеквадратичной ошибки (MSE)

if (mse == 0)
        return 100;

// Расчет PSNR
    int maxPixelValue = originalRange.min+originalRange.max;
    int psnr =  round ( 10.0 * log10((double)(maxPixelValue * maxPixelValue)*count / mse) );

return psnr;
}

template<typename F, typename T>
inline auto process_point(int x, int y, T &&value, F f) -> decltype(f(value),void())
{
   f(value);
}

template<typename F, typename T>
inline auto process_point(int x, int y, T &&value, F f) -> decltype(f(x,y,value), void())
{
   f(x,y,value);
}

template <typename F, typename the_matrix>
auto process_matrix(the_matrix && matrix, F f) -> decltype( matrix[0].size(), matrix.size(), matrix[0][0], void())
{    
    for (int y=0; y<matrix.size();y++)
    for (int x=0; x<matrix[y].size();x++)
        process_point(x,y,matrix[y][x],f);    
}

template <typename F, typename the_matrix>
auto process_matrix(the_matrix &&matrix, Block const& block, F f) -> decltype( matrix[0].size(), matrix.size(), matrix.at(0).at(0), void())
{
    for (int j=0; j < block.y.size; j++)
    for (int i=0; i < block.x.size; i++)
    {
        int x = block.x.offs+i;
        int y = block.y.offs+j;
        process_point(x,y,matrix.at(y).at(x),f);    
    }
}

template <typename F>
inline the_matrix make_matrix(int width, int height, F f)
{
    the_matrix matrix(height, std::vector<int>(width,0));
    process_matrix(matrix, f);
    return matrix;
}

inline the_matrix make_sky(int width, int height) {
    auto result = make_matrix(width, height,
     [sz=width*height](int x, int y, int &v){
        
          if (rand() % 4 == 0)
             v = 255;
         else
            v =  y + x/4;

} );
     cubicBlur3x3(result);
     cubicBlur3x3(result);
     return result;
}

namespace inline_the_matrix {
    
inline int ilog2_32(uint32_t v)
{
   if (v == 0)
      return 1;
   return 32-__builtin_clz(v);
}
}

template <typename the_matrix>
int bitSize(the_matrix && matrix) {
    int bitcount = 0;
    auto range = get_range(matrix);
        
    process_matrix( std::forward<the_matrix>(matrix), [&](int v){        
            
            if ( range.min >= 0)
                bitcount += inline_the_matrix::ilog2_32(v);//unsigned bitsize
            else
                bitcount += inline_the_matrix::ilog2_32(v)+1;//signed bitsize
    });
    
    return bitcount;
}

template <typename the_matrix>
int matrix_energy(the_matrix&& m)
{
    return bitSize(std::forward<the_matrix>(m)) / 8;
}

inline std::vector<int> flatten(const the_matrix & matrix)
{
    std::vector<int> vec; 
    vec.reserve(  matrix.size() * matrix[0].size());
    auto vec_ins = std::back_inserter(vec);
    for (const auto & row : matrix)
    {
        std::copy(row.begin(), row.end(), vec_ins);
    }
    return vec;
}

#endif

Source code

#pragma once
#include <algorithm>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <vector>
#include "utils.hpp"

class Codebook {
    std::vector<int> codebook;
    std::vector<int> value_to_cw;

template <typename Iterator>
    void process_range(Iterator begin, Iterator end, Iterator zero_coefficient) {
        if (begin == end) return;
        auto bin_length = std::distance(begin, end);  // length of bin
        if (bin_length == 1) {
            if (*begin > 0) *begin = -*begin;
            return;
        }

auto maxIt = std::max_element(begin, end);
        *maxIt = -*maxIt;
        if (*maxIt == 0)
            return;
  
     
        auto clear_left = maxIt == begin ? end : maxIt - 1;
        auto clear_right = maxIt + 1;

auto dynamic_bin = std::max<int>(1, std::distance(zero_coefficient, maxIt)/ 4);

for (int i = 1; i < dynamic_bin; i++) {
            if (clear_left != end) {
                //if (std::distance(zero_coefficient,clear_left) >= 4)
                *clear_left = 0;
                // else
                //   *clear_left -= *clear_left;//preserve ditails
                if (clear_left == begin)
                    clear_left = end;
                else
                    --clear_left;
            }
            if (clear_right != end) {
                *clear_right = 0;
                ++clear_right;
            }
        }

if (clear_left != end) process_range(begin, clear_left, zero_coefficient);
        process_range(clear_right, end, zero_coefficient);
    }

void create_codeword_index() {
        value_to_cw.clear();
        for (int i = 0; i <= codebook.back(); i++) {
            value_to_cw.push_back(std::distance(
                codebook.begin(),

find_nerest(codebook.begin(), codebook.end(), i)));
        }
        //std::cerr << "create_codeword_index:" << value_to_cw;
    }

public:
    Codebook(){};
    Codebook(Codebook&& other) noexcept
        : codebook(std::move(other.codebook)),
          value_to_cw(std::move(other.value_to_cw)) {}

Codebook(std::vector<int> hist) {
        //std::cerr << "histogram=" << hist;
     
        // for (auto f4=hist.begin(); f4 != hist.begin() + std::min<int>(5,hist.size()); ++f4)
        // {
        //       if (*f4 > 0)
        //       *f4 = -*f4;
        // }
        process_range(hist.begin(), hist.end(), hist.begin());

for (int i = 0; i < hist.size(); ++i) {
            if (hist[i]) {
                codebook.push_back(i);
            }
        }
        create_codeword_index();
    }
    void set_external_codebook(std::vector<int> external) {
        codebook = std::move(external);
        create_codeword_index();
    }

Codebook(int max_value, int step = 5, int expa = 16) {
        for (auto x = 0, v = 0, progress = 0; x <= max_value && v <= max_value; x++) {
            codebook.push_back(v);
            v += progress / expa + 1;
            progress = progress + step;
        }
        create_codeword_index();
    }

//Codebook() : Codebook(256, 5, 16) {}

int get_codeword(int v) {
        if (is_empty())
            return v; //by-pass
        int max_value = static_cast<int>(value_to_cw.size() - 1);
        if (v <= -max_value) return -value_to_cw.back();
        if (v >= max_value) return value_to_cw.back();
        if (v < 0)
            return -value_to_cw[-v];
        else
            return value_to_cw[v];
    }

int get_value(int cw) {
        if (is_empty())
            return cw; //by-pass
        int max_value = static_cast<int>(codebook.size() - 1);
        if (cw <= -max_value) return -codebook.back();
        if (cw >= max_value) return codebook.back();
        if (cw < 0)
            return -codebook[-cw];
        else
            return codebook[cw];
    }

void clear() {
        codebook.clear();
        value_to_cw.clear();
    }
    bool is_empty() {
        return codebook.size() == 0;
    }
    int size() { return codebook.size(); }
    int at(int cw) { return codebook.at(cw); }
    const std::vector<int>& get_codebook() {
        return codebook;
    }

Codebook& operator=(Codebook&& other) noexcept {
        if (this != &other) {
            codebook = std::move(other.codebook);
            value_to_cw = std::move(other.value_to_cw);
        }
        return *this;
    }
};

Source code

//https://raw.githubusercontent.com/vovach777/huffman_advanced/main/huffman.hpp
#pragma once
#include <algorithm>
#include <bitset>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <limits>
#include <ostream>
#include <queue>
#include <stdexcept>
#include <string>
#include <unordered_map>
#include <vector>
#include "bitstream.hpp"
#include "utils.hpp"

namespace pack {
    namespace impl {
    using namespace std;
    // A Huffman
    template <typename T = char>
    class Huffman {
        // A Huffman node
        struct Node {
            T ch;
            int freq;
            Node *left, *right;
            // Function to allocate a new tree node
            Node(T ch, int freq, Node* left, Node* right)
                : ch(ch), freq(freq), left(left), right(right) {}
            bool operator<(const Node& other) const {
                return freq < other.freq;
            }
        };

// Comparison object to be used to order the heap

struct Node_comp {
            bool operator()(Node* l, Node* r) {
                // highest priority item has lowest frequency
                return l->freq > r->freq;
            }
        };

Node* root;
        vector<Node> pool;
        using BitCountType = uint8_t;
        using HuffmanCodeType = uint32_t;
        unordered_map<T, pair<BitCountType, HuffmanCodeType>> huffmanCode;
        vector<vector<T>> DHT;

// traverse the Huffman Huffman and store Huffman Codes
        // in a map.
        void fillDHT(Node* node, int len) {
            if (node == nullptr) return;

// found a leaf node
            if (!node->left && !node->right) {
                assert(len <= 16 && len > 0);
                if (DHT.size() < len) DHT.resize(len);
                DHT[len - 1].push_back(node->ch);
            }

fillDHT(node->left, len + 1);
            fillDHT(node->right, len + 1);
        }

public:
        void encodeHuffmanTree(BitWriter& dest) {
            if (DHT.size() == 0) throw std::domain_error("No DHT found!");

vector<bool> str;
            for (int i = 0; i < DHT.size(); i++) {
                dest.writeBits(8, DHT[i].size());
                for (auto symbol : DHT[i]) {
                    auto catindex = symbol_to_catindex(symbol);
                    dest.writeBits(4, catindex & 0xf);
                    if ( catindex > 0)
                        dest.writeBits(catindex & 0xf, catindex >> 4);
                }
            }
            // rest of DHT (16 bytes)
            for (int i = DHT.size(); i < 16; i++) {
                dest.writeBits(8, 0);
            }
        }

void decodeHuffmanTree(BitReader& src) {
            // decode huffman table
            DHT.clear();
            DHT.resize(16);

for (int i = 0; i < 16; i++) {
                DHT[i].resize(src.readBits(8));
                if (DHT[i].size() == 0) continue;
                cout << i << ":";
                // cout << "DHT_" << i << ":" << DHT[i].size() << endl;
                for (auto& value : DHT[i]) {
                    auto cat = src.readBits(4);
                    auto index = cat == 0 ? 0 : src.readBits(cat);
                    value = catindex_to_symbol(cat, index);
                    cout << value << ",";
                }
                cout << endl;
            }
            cout << endl;
            create_lockup_table();  // need for encoder
            buildHuffmanTree();     // need for decoder
        }
        void create_lockup_table() {
            //std::cout << "generation huffman codes:" << endl;
            huffmanCode.clear();
            HuffmanCodeType code = 0;
            for (int bits = 0; bits < DHT.size(); ++bits) {
                for (int i = 0; i < DHT[bits].size(); i++) {
                    // HuffmanCodeType rev_code = 0;
                    // for (int n = bits-1; n >= 0; n--){
                    //     rev_code |= (code >> n) & 1;
                    //     rev_code <<= 1;
                    // }

huffmanCode[DHT[bits][i]] = make_pair(bits + 1, code);
                    code++;
                }
                code <<= 1;
            }
            //cout << "Huffman Codes are :\n" << '\n';
            // for (const auto& pair : huffmanCode) {
            //     std::bitset<std::numeric_limits<HuffmanCodeType>::digits> bs(
            //         pair.second.second);
            //     auto s = bs.to_string();
            //     cout << setw(4) << pair.first << " "
            //          << s.substr(s.size() - pair.second.first,
            //                      pair.second.first)
            //          << " (" << (int)pair.second.first << ","
            //          << pair.second.second << ")" << endl;
            // }
        }
        // Builds Huffman Huffman
        template <typename Iterator>
        auto buildHuffmanTree(Iterator begin, Iterator end)
            -> decltype(std::enable_if_t<
                            std::is_same_v<typename std::iterator_traits<
                                               Iterator>::value_type,
                                           T>,
                            Iterator>(),
                        void()) {
            // count frequency of appearance of each character
            // and store it in a map
            unordered_map<T, int> freq;
            for (auto it = begin; it != end; ++it) {
                freq[*it]++;
            }

pool.clear();
            pool.reserve(freq.size() * 2);

// Create a priority queue to store live nodes of
            // Huffman tree;

priority_queue<Node*, vector<Node*>, Node_comp> pq;

// Create a leaf node for each character and add it
            // to the priority queue.
            for (auto pair : freq) {
                pq.push(&pool.emplace_back(pair.first, pair.second, nullptr,
                                           nullptr));
            }

// do till there is more than one node in the queue
            while (pq.size() != 1) {
                // Remove the two nodes of highest priority
                // (lowest frequency) from the queue
                Node* left = pq.top();
                pq.pop();
                Node* right = pq.top();
                pq.pop();

// Create a new internal node with these two nodes
                // as children and with frequency equal to the sum
                // of the two nodes' frequencies. Add the new node
                // to the priority queue.
                int sum = left->freq + right->freq;
                pq.push(&pool.emplace_back('\0', sum, left, right));
            }

// root stores pointer to root of Huffman Huffman
            root = pq.top();

// traverse the Huffman Huffman and store (not Huffman Codes) bitlength
            // of each symbol in a map. Also prints them
            DHT.clear();
            DHT.reserve(16);
            fillDHT(root, 0);

// huffman table not needed any more - lockup table to be
            // regenerated from DHT
            root = nullptr;
            pool.clear();
            create_lockup_table();
            buildHuffmanTree();
        }
        // from lockuptable
        void buildHuffmanTree() {
            if (huffmanCode.size() == 0) {
                root = nullptr;
                return;
            }
            pool.clear();
            pool.reserve(huffmanCode.size() * 2);
            root = &pool.emplace_back('\0', 1, nullptr, nullptr);

// for (const auto & symbols : DHT)
            // for (const auto symbol : symbols) {
            //     const auto& code = huffmanCode[symbol];
            //     cout << symbol << " " << code<< '\n';
            for (const auto& pair : huffmanCode) {
                auto symbol = pair.first;
                // cerr << endl << "(" << symbol << ") : ";
                auto& code = pair.second;
                auto n = root;
                for (int i = code.first - 1; i >= 0; --i) {
                    auto v = (code.second >> i) & 1;
                    // cerr << v;
                    if (v == 0) {
                        if (!n->left) {
                            n->left =
                                &pool.emplace_back('\0', 1, nullptr, nullptr);
                        }
                        assert(n->freq == 1);
                        n = n->left;

} else {
                        if (!n->right) {
                            n->right =
                                &pool.emplace_back('\0', 1, nullptr, nullptr);
                        }
                        assert(n->freq == 1);
                        n = n->right;
                    }
                }
                assert(n->left == nullptr && n->right == nullptr &&
                       n->freq == 1);
                n->ch = symbol;
                n->freq = 0;
            }
            // cerr << endl;
        }

template <typename Iterator>
        void encode(BitWriter& dest, Iterator begin, Iterator end) {
            for (auto it = begin; it != end; ++it) {
                encode(dest, *it);
            }
        }
        inline void encode(BitWriter& dest, T value) {
            auto [bit_count, code] = huffmanCode.at(value);
            // while (bit_count--) {
            //     dest.writeBit(code & 1);
            //     code >>= 1;
            // }
            dest.writeBits(bit_count, code);
        }

// traverse the Huffman Huffman and decode the encoded string
        T decode(BitReader& src) {
            if (root == nullptr) {
                throw std::domain_error("No the huffman-tree found!");
            }

auto n = root;

for (;;) {
                // if (src.bit_left() == 0)
                //     return root->ch;

if (src.readBit() == false) {
                    n = n->left;
                } else {
                    n = n->right;
                }
                // found a leaf node
                if (n->freq == 0) {
                    // cerr << "(" << n->ch << ")";
                    return n->ch;
                }
            }
        }
        template <typename Iterator>
        void decode(BitReader& src, Iterator begin, Iterator end) {
            for (auto it = begin; it < end; it++) {
                *it = decode(src);
            }
        }
    };
    }
    //export section
    template <typename T = char>
    using Huffman = impl::Huffman<T>;
}