diff --git a/CMakeLists.txt b/CMakeLists.txt index 71c5e7b..a4f9297 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,14 +1,16 @@ cmake_minimum_required (VERSION 2.6) +project(compress) include_directories(SYSTEM "${PROJECT_SOURCE_DIR}/include/") -add_definitions( -std=c++0x ) +set(CMAKE_CXX_FLAGS ${CMAKE_CXX_FLAGS} "-std=c++11") -option(MultiCoreBSP "Building with Multi Core BSP" OFF) +option(MultiCoreBSP "Building with Multi Core BSP" ON) if(MultiCoreBSP) # setup for my mac set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS}") - set(libs mcbsp1.1.0 m pthread) + set(libs mcbsp m pthread) add_definitions( -DUSE_MCBSP ) + include_directories(SYSTEM "${PROJECT_SOURCE_DIR}/contrib/mcbsp") else() # setup for cartesius set(libs bsponmpi m) @@ -19,3 +21,4 @@ include_directories(SYSTEM ${Boost_INCLUDE_DIRS}) set(libs ${libs} ${Boost_LIBRARIES}) add_subdirectory("wavelet") +add_subdirectory("contrib/mcbsp") diff --git a/wavelet/jcmp.cpp b/wavelet/jcmp.cpp index 183c772..4305995 100644 --- a/wavelet/jcmp.cpp +++ b/wavelet/jcmp.cpp @@ -1,21 +1,38 @@ #include #include +#include "jcmp.hpp" +#include "wavelet.hpp" +#include "remap.hpp" + #include +#include #include -#include "jcmp.hpp" -#include "wavelet.hpp" +#include +#include using namespace wvlt; using namespace jcmp; // note: we take a copy, because we will modify it in place -static image compress(std::vector img, uint width, double threshold, unsigned int& nzeros){ - uint height = img.size() / width; +static image compress(std::vector const & img0, uint16_t width, double threshold, unsigned int& nzeros){ + uint16_t height = img0.size() / width; assert(is_pow_of_two(width)); assert(is_pow_of_two(height)); + assert(width == height); + + std::vector img(img0.size(), 0); + for(uint16_t y = 0; y < height; y++){ + for(uint16_t x = 0; x < width; x++) { + img[remap::to_hilbert(width, x, y)] = img0[x + y*width]; + } + } + + wavelet(&img[0], width * height, 1); + +#if 0 // wavelet transform in x-direction for(unsigned int i = 0; i < height; ++i){ wavelet(&img[i*width], width, 1); @@ -25,6 +42,7 @@ static image compress(std::vector img, uint width, double threshold, uns for(unsigned int i = 0; i < width; ++i){ wavelet(&img[i], height, width); } +#endif double min_abs = 10000.0; double max_abs = 0.0; @@ -43,8 +61,8 @@ static image compress(std::vector img, uint width, double threshold, uns // save the principal coefficients std::vector v; - for(unsigned int y = 0; y < height; ++y){ - for(unsigned int x = 0; x < width; ++x){ + for(uint16_t y = 0; y < height; ++y){ + for(uint16_t x = 0; x < width; ++x){ auto&& el = img[x + width*y]; if(el != 0) v.push_back({q.forwards(el), uint(x), uint(y)}); } @@ -55,8 +73,8 @@ static image compress(std::vector img, uint width, double threshold, uns } static std::vector decompress(image in){ - auto width = in.header.width; - auto height = in.header.height; + const auto width = in.header.width; + const auto height = in.header.height; assert(is_pow_of_two(width)); assert(is_pow_of_two(height)); @@ -66,12 +84,13 @@ static std::vector decompress(image in){ // read in coefficient on coordinates for(auto it = in.data.begin(); it != in.data.end(); ++it){ - auto&& x = *it; - image[x.x + width*x.y] = q.backwards(x.c); + auto&& coef = *it; + image[coef.x + width*coef.y] = q.backwards(coef.c); } in.clear(); +#if 0 // inverse wavelet transform in y-direction for(unsigned int i = 0; i < width; ++i){ unwavelet(&image[i], height, width); @@ -81,18 +100,70 @@ static std::vector decompress(image in){ for(unsigned int i = 0; i < height; ++i){ unwavelet(&image[i*width], width, 1); } +#endif + + unwavelet(&image[0], width*height, 1); + + std::vector image2(width * height, 0.0); + + for(uint32_t d = 0; d < width * height; ++d) { + const auto p = remap::from_hilbert(width, d); + image2[p.first + width * p.second] = image[d]; + } - return image; + return image2; } -int main(){ +int main(int argc, char* argv[]){ + using namespace std; namespace fs = boost::filesystem; - fs::path directory("images"); - fs::directory_iterator eod; - for(fs::directory_iterator it(directory); it != eod; ++it){ - auto && path = it->path(); - if(path.extension() != ".png") continue; + struct { + vector filenames = {}; + double threshold = 5.0; + } options; + + namespace po = boost::program_options; + // Describe program options + po::options_description opts; + opts.add_options() + ("threshold", po::value(), "") + ("help", po::bool_switch(), "show this help") + ("file", po::value< vector >(), "input file(s)"); + // and inputs + po::positional_options_description p; + p.add("file", -1); + + // Parse and set options + try { + po::variables_map vm; + po::store(po::command_line_parser(argc, argv). + options(opts).positional(p).run(), vm); + po::notify(vm); + + if(vm["help"].as()){ + std::cout << "jcmp" << std::endl; + std::cout << opts << std::endl; + return 0; + } + + if(vm.count("file")){ + options.filenames = vm["file"].as< vector >(); + } + + if(vm.count("threshold")) { + options.threshold = vm["threshold"].as(); + } + } catch(std::exception& e){ + std::cout << colors::red("ERROR: ") << e.what() << std::endl; + std::cout << opts << std::endl; + return 1; + } + + for(const auto & path : options.filenames) { + if(path.extension() != ".png") { + cerr << "WARNING: can only read .png, but " << path << " is not" << endl; + } // open file std::string filename = path.string(); @@ -109,7 +180,7 @@ int main(){ // compress and decompress to see how we lost information unsigned int nzeros = 0; - auto compressed_vec = decompress(compress(image_vec, width, 5.0, nzeros)); + auto compressed_vec = decompress(compress(image_vec, width, options.threshold, nzeros)); // output some information std::cout << field("raw") << width*height << std::endl; @@ -117,7 +188,7 @@ int main(){ std::cout << field("nz") << nzeros << std::endl; // save to output file - std::string cfilename = "compressed/" + std::to_string(nzeros) + path.filename().string(); + std::string cfilename = "hil_compressed_" + std::to_string(nzeros) + path.filename().string(); png::gray_ostream compressed_image(width, height, cfilename); for(unsigned int i = 0; i < compressed_vec.size(); ++i) compressed_image << compressed_vec[i]; } diff --git a/wavelet/jcmp_quantization.hpp b/wavelet/jcmp_quantization.hpp index b904edc..fb76fc1 100644 --- a/wavelet/jcmp_quantization.hpp +++ b/wavelet/jcmp_quantization.hpp @@ -50,6 +50,7 @@ namespace jcmp { enum class type { linear, logarithmic, + quadric_root, square_root }; @@ -71,6 +72,10 @@ namespace jcmp { b.f = &log; b.f_inv = &exp; break; + case type::quadric_root: + b.f = [](double x){ return 1 - 1/(0.25*x*x + 1); }; + b.f_inv = [](double x) { return std::sqrt(1/(1 - x)-1); }; + break; case type::square_root: b.f = &sqrt; b.f_inv = &sqr; diff --git a/wavelet/remap.hpp b/wavelet/remap.hpp new file mode 100644 index 0000000..4ab9545 --- /dev/null +++ b/wavelet/remap.hpp @@ -0,0 +1,100 @@ +#pragma once + +#include +#include +#include + +namespace remap { +// there are some more efficient ways in +// http://graphics.stanford.edu/~seander/bithacks.html#InterleaveTableObvious +// but I was too lazy to understand them + +// some templates dealing with sizes +// clang-format off +template struct unsigned_size {}; +template <> struct unsigned_size { using twice = uint16_t; /* 4 bits */ }; +template <> struct unsigned_size { using twice = uint32_t; using half = uint8_t; }; +template <> struct unsigned_size { using twice = uint64_t; using half = uint16_t; }; +template <> struct unsigned_size { /* 128 bits */ using half = uint32_t; }; +// clang-format on + +// shortcuts +template +using twice = typename unsigned_size::twice; +template +using half = typename unsigned_size::half; + +// converts x,y coordinates to z order +template +twice to_z_order(const T xin, const T yin) { + const twice x = xin; + const twice y = yin; + twice z = 0; + for (int i = 0; i < sizeof(T) * CHAR_BIT; i++) { + z |= (x & (1U << i)) << i | (y & (1U << i)) << (i + 1); + } + return z; +} + +// converts z to x,y +template +std::pair, half> from_z_order(const T z) { + half x = 0; + half y = 0; + for (int i = 0; i < sizeof(T) * CHAR_BIT; i += 2) { + x |= (z & (1U << i)) >> (i / 2); + y |= (z & (1U << (i + 1))) >> (i / 2 + 1); + } + return {x, y}; +} + + +// from wikipedia +// rotate/flip a quadrant appropriately +template +void hilbert_rot(T n, T * x, T * y, T rx, T ry) { + if (ry == 0) { + if (rx == 1) { + *x = n - 1 - *x; + *y = n - 1 - *y; + } + + // Swap x and y + int t = *x; + *x = *y; + *y = t; + } +} + +template +twice to_hilbert(twice n, T xin, T yin) { + twice x = xin; + twice y = yin; + twice d = 0; + for (twice s = n / 2; s > 0; s /= 2) { + twice rx = (x & s) > 0; + twice ry = (y & s) > 0; + d += s * s * ((3 * rx) ^ ry); + hilbert_rot(s, &x, &y, rx, ry); + } + return d; +} + +// convert d to (x,y) +template +std::pair, half> from_hilbert(T n, T d) { + T t = d; + T x = 0; + T y = 0; + for (T s = 1; s < n; s *= 2) { + T rx = 1 & (t / 2); + T ry = 1 & (t ^ rx); + hilbert_rot(s, &x, &y, rx, ry); + x += s * rx; + y += s * ry; + t /= 4; + } + return {half(x), half(y)}; +} + +} diff --git a/wavelet/remap_test.cpp b/wavelet/remap_test.cpp new file mode 100644 index 0000000..152b44c --- /dev/null +++ b/wavelet/remap_test.cpp @@ -0,0 +1,66 @@ +#include "remap.hpp" + +#include +#include + +using namespace std; + +template +static void test_z_order(uint64_t xmax, uint64_t ymax){ + uint64_t c = 10; + for (uint64_t x2 = 0; x2 <= xmax; ++x2) { + for (uint64_t y2 = 0; y2 <= ymax; ++y2) { + const T x = x2; + const T y = y2; + const auto z = remap::to_z_order(x, y); + const auto p = remap::from_z_order(z); + + if (x != p.first || y != p.second) { + cout << +x << ' ' << +y << " -> " << z << " -> " << +p.first << ' ' << +p.second << '\n'; + c--; + if(!c) return; + } + } + } +} + +template +static void test_z_order_max(){ + return test_z_order(std::numeric_limits::max(), std::numeric_limits::max()); +} + +template +static void test_hilbert(uint64_t xmax){ + uint64_t c = 10; + for (uint64_t x2 = 0; x2 <= xmax; ++x2) { + for (uint64_t y2 = 0; y2 <= xmax; ++y2) { + const T x = x2; + const T y = y2; + const remap::twice xmax2 = xmax+1; + const auto z = remap::to_hilbert(xmax2, x, y); + const auto p = remap::from_hilbert(xmax2, z); + + if (x != p.first || y != p.second) { + cout << +x << ' ' << +y << " -> " << z << " -> " << +p.first << ' ' << +p.second << '\n'; + c--; + if(!c) return; + } + } + } +} + +template +static void test_hilbert_max(){ + return test_hilbert(std::numeric_limits::max()); +} + +int main(int argc, char * argv[]) { + //cout << "Z order: 8" << endl; + //test_z_order_max(); + cout << "Hilbert: 8" << endl; + test_hilbert_max(); + //cout << "Z order: 16" << endl; + //test_z_order_max(); + cout << "Hilbert: 16" << endl; + test_hilbert_max(); +}