Adds hilbert/wavelet thingies.

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")

wavelet/jcmp.cpp

@@ -1,21 +1,38 @@
 #include <includes.hpp>
 #include <utilities.hpp>
 
+#include "jcmp.hpp"
+#include "wavelet.hpp"
+#include "remap.hpp"
+
 #include <boost/filesystem.hpp>
+#include <boost/program_options.hpp>
 #include <png.hpp>
 
-#include "jcmp.hpp"
+#include <vector>
-#include "wavelet.hpp"
+#include <iostream>
 
 using namespace wvlt;
 using namespace jcmp;
 
 // note: we take a copy, because we will modify it in place
-static image compress(std::vector<double> img, uint width, double threshold, unsigned int& nzeros){
+static image compress(std::vector<double> const & img0, uint16_t width, double threshold, unsigned int& nzeros){
-	uint height = img.size() / width;
+	uint16_t height = img0.size() / width;
 	assert(is_pow_of_two(width));
 	assert(is_pow_of_two(height));
+	assert(width == height);
+
+	std::vector<double> 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<double> 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<double> img, uint width, double threshold, uns
 
 	// save the principal coefficients
 	std::vector<jcmp::coefficient> v;
-	for(unsigned int y = 0; y < height; ++y){
+	for(uint16_t y = 0; y < height; ++y){
-		for(unsigned int x = 0; x < width; ++x){
+		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<double> img, uint width, double threshold, uns
 }
 
 static std::vector<double> decompress(image in){
-	auto width = in.header.width;
+	const auto width = in.header.width;
-	auto height = in.header.height;
+	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<double> decompress(image in){
 
 	// read in coefficient on coordinates
 	for(auto it = in.data.begin(); it != in.data.end(); ++it){
-		auto&& x = *it;
+		auto&& coef = *it;
-		image[x.x + width*x.y] = q.backwards(x.c);
+		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<double> 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<double> 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");
+	struct {
-	fs::directory_iterator eod;
+		vector<fs::path> filenames = {};
-	for(fs::directory_iterator it(directory); it != eod; ++it){
+		double threshold = 5.0;
-		auto && path = it->path();
+	} options;
-		if(path.extension() != ".png") continue;
+
+	namespace po = boost::program_options;
+	// Describe program options
+	po::options_description opts;
+	opts.add_options()
+		("threshold", po::value<double>(), "")
+		("help", po::bool_switch(), "show this help")
+		("file", po::value< vector<fs::path> >(), "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<bool>()){
+			std::cout << "jcmp" << std::endl;
+			std::cout << opts << std::endl;
+			return 0;
+		}
+
+		if(vm.count("file")){
+			options.filenames = vm["file"].as< vector<fs::path> >();
+		}
+
+		if(vm.count("threshold")) {
+			options.threshold = vm["threshold"].as<double>();
+		}
+	} 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];
 	}

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;

wavelet/remap.hpp

@@ -0,0 +1,100 @@
+#pragma once
+
+#include <climits>
+#include <cstdint>
+#include <utility>
+
+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 <typename T> struct unsigned_size {};
+template <> struct unsigned_size<uint8_t>  { using twice = uint16_t;      /* 4 bits */      };
+template <> struct unsigned_size<uint16_t> { using twice = uint32_t; using half = uint8_t;  };
+template <> struct unsigned_size<uint32_t> { using twice = uint64_t; using half = uint16_t; };
+template <> struct unsigned_size<uint64_t> {      /* 128 bits */     using half = uint32_t; };
+// clang-format on
+
+// shortcuts
+template <typename T>
+using twice = typename unsigned_size<T>::twice;
+template <typename T>
+using half = typename unsigned_size<T>::half;
+
+// converts x,y coordinates to z order
+template <typename T>
+twice<T> to_z_order(const T xin, const T yin) {
+	const twice<T> x = xin;
+	const twice<T> y = yin;
+	twice<T> 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 <typename T>
+std::pair<half<T>, half<T>> from_z_order(const T z) {
+	half<T> x = 0;
+	half<T> 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 <typename T>
+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 <typename T>
+twice<T> to_hilbert(twice<T> n, T xin, T yin) {
+	twice<T> x = xin;
+	twice<T> y = yin;
+	twice<T> d = 0;
+	for (twice<T> s = n / 2; s > 0; s /= 2) {
+		twice<T> rx = (x & s) > 0;
+		twice<T> 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 <typename T>
+std::pair<half<T>, half<T>> 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<T>(x), half<T>(y)};
+}
+
+}

wavelet/remap_test.cpp

@@ -0,0 +1,66 @@
+#include "remap.hpp"
+
+#include <iostream>
+#include <limits>
+
+using namespace std;
+
+template <typename T>
+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 <typename T>
+static void test_z_order_max(){
+	return test_z_order<T>(std::numeric_limits<T>::max(), std::numeric_limits<T>::max());
+}
+
+template <typename T>
+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<T> 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 <typename T>
+static void test_hilbert_max(){
+	return test_hilbert<T>(std::numeric_limits<T>::max());
+}
+
+int main(int argc, char * argv[]) {
+	//cout << "Z order: 8" << endl;
+	//test_z_order_max<uint8_t>();
+	cout << "Hilbert: 8" << endl;
+	test_hilbert_max<uint8_t>();
+	//cout << "Z order: 16" << endl;
+	//test_z_order_max<uint16_t>();
+	cout << "Hilbert: 16" << endl;
+	test_hilbert_max<uint16_t>();
+}