diff --git a/wavelet/jcmp.cpp b/wavelet/jcmp.cpp index 1fde5e7..b768350 100644 --- a/wavelet/jcmp.cpp +++ b/wavelet/jcmp.cpp @@ -7,40 +7,29 @@ #include "jcmp.hpp" #include "wavelet.hpp" -using namespace wvlt::V2; - -// Can be any functions from R+ -> R -// as long as backward is its inverse -static double forward(double x){ - //return std::log(1 + std::log(1 + x)); - return std::log(x); -} - -static double backward(double x){ - //return std::exp(std::exp(x) - 1) - 1; - return std::exp(x); -} +using namespace wvlt; +using namespace jcmp; // note: we take a copy, because we will modify it in place -static jcmp::image compress(std::vector image, jcmp::uint width, double threshold, unsigned int& nzeros){ - jcmp::uint height = image.size() / width; +static image compress(std::vector img, uint width, double threshold, unsigned int& nzeros){ + uint height = img.size() / width; assert(is_pow_of_two(width)); assert(is_pow_of_two(height)); // wavelet transform in x-direction for(unsigned int i = 0; i < height; ++i){ - wavelet(&image[i*width], width, 1); + wavelet(&img[i*width], width, 1); } // wavelet transform in y-direction for(unsigned int i = 0; i < width; ++i){ - wavelet(&image[i], height, width); + wavelet(&img[i], height, width); } double min_abs = 10000.0; double max_abs = 0.0; - for(auto& el : image){ + for(auto& el : img){ auto absel = std::abs(el); if(absel > threshold) { min_abs = std::min(min_abs, absel); @@ -50,28 +39,28 @@ static jcmp::image compress(std::vector image, jcmp::uint width, double } } - jcmp::quantization q(&forward, &backward, max_abs, min_abs); + auto q = quantization::get(quantization::type::logarithmic, max_abs, min_abs); // save the principal coefficients std::vector v; for(unsigned int y = 0; y < height; ++y){ for(unsigned int x = 0; x < width; ++x){ - auto&& el = image[x + width*y]; - if(el != 0) v.push_back({q.forwards(el), jcmp::uint(x), jcmp::uint(y)}); + auto&& el = img[x + width*y]; + if(el != 0) v.push_back({q.forwards(el), uint(x), uint(y)}); } } nzeros = v.size(); - return jcmp::image(width, height, q.p, std::move(v)); + return image(width, height, q.p, std::move(v)); } -static std::vector decompress(jcmp::image in){ +static std::vector decompress(image in){ auto width = in.header.width; auto height = in.header.height; assert(is_pow_of_two(width)); assert(is_pow_of_two(height)); - auto q = in.header.get_quantization(&forward, &backward); + auto q = in.header.get_quantization(quantization::type::logarithmic); std::vector image(width * height, 0.0); @@ -120,7 +109,7 @@ int main(){ // compress and decompress to see how we lost information unsigned int nzeros = 0; - auto compressed_vec = decompress(compress(image_vec, width, 0.1, nzeros)); + auto compressed_vec = decompress(compress(image_vec, width, 0.5, nzeros)); // output some information std::cout << field("raw") << human_string(sizeof(uint8_t) * image_vec.size(), "b") << std::endl; diff --git a/wavelet/jcmp_image.hpp b/wavelet/jcmp_image.hpp index 0094412..8c2ecdd 100644 --- a/wavelet/jcmp_image.hpp +++ b/wavelet/jcmp_image.hpp @@ -13,10 +13,10 @@ namespace jcmp { uint width; uint height; uint length; - quantize_params qp; + quantization::parameters qp; header() = default; - header(uint width_, uint height_, uint length_, quantize_params const & p) + header(uint width_, uint height_, uint length_, quantization::parameters const & p) : signature{'J', 'C', 'M', 'P'} , width(width_) , height(height_) @@ -24,13 +24,8 @@ namespace jcmp { , qp(p) {} - template - quantization get_quantization(F const & f, F const & f_inv){ - quantization q; - q.p = qp; - q.f = f; - q.f_inv = f_inv; - return q; + quantization::base get_quantization(quantization::type t){ + return quantization::get(t, qp.f_max_abs, qp.f_min_abs, false); } }; @@ -46,7 +41,7 @@ namespace jcmp { image() = default; - image(uint width, uint height, quantize_params const & p, std::vector const & data_in) + image(uint width, uint height, quantization::parameters const & p, std::vector const & data_in) : header(width, height, data_in.size(), p) , data(data_in) {} diff --git a/wavelet/jcmp_quantization.hpp b/wavelet/jcmp_quantization.hpp index fbf2a4c..3f06ca2 100644 --- a/wavelet/jcmp_quantization.hpp +++ b/wavelet/jcmp_quantization.hpp @@ -1,49 +1,91 @@ #pragma once -namespace jcmp { - // input in [0, 1], output [0, 255] - inline uint8_t clamp_round(double x){ - if(x <= 0.0) return 0; - if(x >= 1.0) return 255; - return uint8_t(std::round(255.0*x)); - } +#include - // parameters which need to be stored in the file - // to be able to dequantize - struct quantize_params { - double f_max_abs; - double f_min_abs; - }; - - // Quantization may be non linear (given by f and f_inv) - struct quantization { - std::function f; - std::function f_inv; - quantize_params p; - - quantization() = default; - - template - quantization(F const & f_, F const & f_inv_, double max_abs, double min_abs) - : f(f_) - , f_inv(f_inv_) - , p{f(max_abs), f(min_abs)} - {} - - uint8_t forwards(double x){ - double y = std::abs(x); - y = (f(y) - p.f_min_abs) / (p.f_max_abs - p.f_min_abs); - if(x < 0) y = -y; - return clamp_round(0.5 * (y + 1.0)); +namespace jcmp { + // Will quantize a double such that it fits in a uint8_t + // Also support non-linear quantization (which is usefull + // when the data is non-uniformly distributed) + namespace quantization { + // input in [0, 1], output [0, 255] + inline uint8_t clamp_round(double x){ + if(x <= 0.0) return 0; + if(x >= 1.0) return 255; + return uint8_t(std::round(255.0*x)); } - double backwards(uint8_t x){ - double y = 2.0 * x / 255.0 - 1.0; - y = (p.f_max_abs - p.f_min_abs) * y; - if(y < 0) return -f_inv(-y + p.f_min_abs); - return f_inv(y + p.f_min_abs); + // parameters which need to be stored in the file + // to be able to dequantize + // TODO: add type of quantization + struct parameters { + double f_max_abs; + double f_min_abs; + }; + + // the basic struct doing the quantization + struct base { + // f needs to be of type R+ -> R + // and F-inv should be its inverse + std::function f; + std::function f_inv; + parameters p; + + uint8_t forwards(double x){ + double y = std::abs(x); + y = (f(y) - p.f_min_abs) / (p.f_max_abs - p.f_min_abs); + if(x < 0) y = -y; + return clamp_round(0.5 * (y + 1.0)); + } + + double backwards(uint8_t x){ + double y = 2.0 * x / 255.0 - 1.0; + y = (p.f_max_abs - p.f_min_abs) * y; + if(y < 0) return -f_inv(-y + p.f_min_abs); + return f_inv(y + p.f_min_abs); + } + }; + + // kinds of quantization which come out of the + // box with get(...) + enum class type { + linear, + logarithmic, + square_root + }; + + // non-overloaded versions + double id(double x) { return x; } + double log(double x) { return std::log(x); } + double exp(double x) { return std::exp(x); } + double sqrt(double x) { return std::sqrt(x); } + double sqr(double x) { return x*x; } + + base get(type t, double max_abs, double min_abs, bool apply = true){ + base b; + switch (t) { + case type::linear: + b.f = &id; + b.f_inv = &id; + break; + case type::logarithmic: + b.f = &log; + b.f_inv = &exp; + break; + case type::square_root: + b.f = &sqrt; + b.f_inv = &sqr; + break; + } + if(apply){ + b.p.f_max_abs = b.f(max_abs); + b.p.f_min_abs = b.f(min_abs); + } else { + b.p.f_max_abs = max_abs; + b.p.f_min_abs = min_abs; + } + return b; } - }; - static_assert(sizeof(quantize_params) == 16, "struct not propery packed"); + static_assert(sizeof(parameters) == 16, "struct not propery packed"); + } }