Last active
December 17, 2015 16:19
-
-
Save YaLTeR/174dd504843d7e5c6268 to your computer and use it in GitHub Desktop.
Implementation of the fractal image compression algorithm. Can compress using fixed block size (4x4, 8x8, 16x16) as well as using a quadtree with a quality setting (1 to 100 inclusive). This is not the complete source code of a program, however the algorithm is here completely.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #include "stdafx.h" | |
| #include <algorithm> | |
| #include <cmath> | |
| #include <cstdint> | |
| #include <cstdio> | |
| #include <exception> | |
| #include <utility> | |
| #include "colorspaces.hpp" | |
| #include "psnr.hpp" | |
| #include "compress\compress.h" | |
| #include "bitarr.hpp" | |
| enum transformation | |
| { | |
| NONE = 0, | |
| ROT_90, | |
| ROT_180, | |
| ROT_270, | |
| FLIP, | |
| FLIP_90, | |
| FLIP_180, | |
| FLIP_270 | |
| }; | |
| struct square | |
| { | |
| std::vector<double> data; | |
| square flipped() const | |
| { | |
| size_t s = static_cast<size_t>(sqrt(data.size())); | |
| square res; | |
| res.data.resize(data.size()); | |
| for (size_t i = 0; i < s; ++i) { | |
| for (size_t j = 0; j < s; ++j) { | |
| res.data[i * s + j] = data[i * s + s - j - 1]; | |
| } | |
| } | |
| return res; | |
| } | |
| square transformed(transformation t) const | |
| { | |
| size_t s = static_cast<size_t>(sqrt(data.size())); | |
| auto tr = static_cast<unsigned char>(t); | |
| square temp; | |
| if (tr >= FLIP) { | |
| temp = flipped(); | |
| tr -= FLIP; | |
| } else { | |
| temp = *this; | |
| } | |
| if (tr == NONE) | |
| return temp; | |
| square res; | |
| res.data.resize(data.size()); | |
| for (size_t i = 0; i < s; ++i) { | |
| for (size_t j = 0; j < s; ++j) { | |
| if (tr == ROT_90) | |
| res.data[j * s + s - i - 1] = temp.data[i * s + j]; | |
| else if (tr == ROT_180) | |
| res.data[(s - i - 1) * s + s - j - 1] = temp.data[i * s + j]; | |
| else if (tr == ROT_270) | |
| res.data[(s - j - 1) * s + i] = temp.data[i * s + j]; | |
| } | |
| } | |
| return res; | |
| } | |
| }; | |
| struct mapping | |
| { | |
| mapping(int idx_b, transformation tr, double q_, double p_, double dist_ = 0.0) | |
| : idx_big(idx_b) | |
| , t(tr) | |
| , q(q_) | |
| , p(p_) | |
| , dist(dist_) | |
| {} | |
| int idx_big; | |
| transformation t; | |
| double q; | |
| double p; | |
| double dist; | |
| }; | |
| struct quadtree_node | |
| { | |
| quadtree_node(size_t bs) | |
| : my_block_size(bs) | |
| {} | |
| bool leftup_is_quadtree; | |
| bool rightup_is_quadtree; | |
| bool leftdown_is_quadtree; | |
| bool rightdown_is_quadtree; | |
| size_t my_block_size; | |
| int x; | |
| int y; | |
| union { | |
| quadtree_node *node; | |
| mapping *map; | |
| } lu; | |
| union { | |
| quadtree_node *node; | |
| mapping *map; | |
| } ru; | |
| union { | |
| quadtree_node *node; | |
| mapping *map; | |
| } ld; | |
| union { | |
| quadtree_node *node; | |
| mapping *map; | |
| } rd; | |
| }; | |
| struct mapping_parent | |
| { | |
| quadtree_node *parent; | |
| enum pos { | |
| LU, RU, LD, RD | |
| } position; | |
| mapping_parent(quadtree_node* par_, pos pos_) | |
| : parent(par_) | |
| , position(pos_) | |
| {} | |
| }; | |
| std::vector<square> partition(const std::vector<double>& a, int w, int block_size) | |
| { | |
| assert(w % block_size == 0); | |
| int w_res = w / block_size; | |
| std::vector<square> res(w_res * w_res); | |
| for (int y = 0; y < w; ++y) { | |
| for (int x = 0; x < w; ++x) { | |
| res[(y / block_size) * w_res + x / block_size].data.push_back(a[y * w + x]); | |
| } | |
| } | |
| return res; | |
| } | |
| void add_transforms(std::vector<square>& a) | |
| { | |
| size_t s = a.size(); | |
| a.reserve(a.size() * 8); | |
| for (int i = 1; i < 8; ++i) { | |
| for (size_t j = 0; j < s; ++j) { | |
| a.push_back(a[j].transformed(static_cast<transformation>(i))); | |
| } | |
| } | |
| } | |
| static const constexpr auto Q_MAX = 0.9, Q_MIN = -1.0; | |
| void compute_best_qp(const square& small, const square& big, size_t block_size, size_t block_size_big_mul, double& q, double& p) | |
| { | |
| size_t s_small = block_size; | |
| size_t s_big = block_size * block_size_big_mul; | |
| double dij = 0.0, dijsq = 0.0, rij = 0.0, rijdij = 0.0; | |
| for (size_t j = 0; j < s_small; ++j) { | |
| for (size_t i = 0; i < s_small; ++i) { | |
| auto idx = j * s_small + i; | |
| auto temp = big.data[j * block_size_big_mul * s_big + i * block_size_big_mul]; | |
| dij += temp; | |
| dijsq += temp * temp; | |
| rij += small.data[idx]; | |
| rijdij += small.data[idx] * temp; | |
| } | |
| } | |
| q = (rijdij - dij * rij) / (dijsq - dij * dij); | |
| //q = 0.75; | |
| q = (q > Q_MAX) ? Q_MAX : q; | |
| q = (q < Q_MIN) ? Q_MIN : q; | |
| // Account for the size reduction we're gonna do later. | |
| auto temp = (q - Q_MIN) * (255 / (Q_MAX - Q_MIN)); | |
| if (temp > 255.0) | |
| temp = 255.0; | |
| else if (temp < 0.0) | |
| temp = 0.0; | |
| auto temp2 = static_cast<uint8_t>(temp); | |
| q = temp2 * ((Q_MAX - Q_MIN) / 255) + Q_MIN; | |
| // Now compute p. | |
| p = (rij - q * dij) / (block_size * block_size); | |
| //p = (rij - dij * (rijdij / dijsq)) / (1 - dij * (dij / dijsq)); | |
| //p = (p <= -1) ? -1 : p; | |
| //p = (p >= 1) ? 1 : p; | |
| //q = (rijdij - p * dij) / dijsq; | |
| //printf("q = %.8f; p = %.8f\n", q, p); | |
| } | |
| std::vector<mapping> make_mapping(const std::vector<square>& par_small, const std::vector<square>& par_big, size_t block_size, size_t block_size_big_mul) | |
| { | |
| int s = par_small.size(); | |
| size_t s_b = par_big.size() / 8; | |
| std::vector<mapping> res(s, mapping(0, NONE, 0.0, 0.0)); | |
| #pragma omp parallel for | |
| for (int i = 0; i < s; ++i) { | |
| //printf("%d\n", i); | |
| auto sq_small = par_small[i]; | |
| // Find best mapping. | |
| int best_idx; | |
| transformation best_t; | |
| double best_q, best_p; | |
| double best_dist = DBL_MAX; | |
| for (size_t j = 0; j < par_big.size(); ++j) { | |
| auto sq_big = par_big[j]; | |
| double q, p; | |
| compute_best_qp(sq_small, sq_big, block_size, block_size_big_mul, q, p); | |
| double dist = 0.0; | |
| for (size_t a = 0; a < block_size; ++a) { | |
| for (size_t b = 0; b < block_size; ++b) { | |
| auto t = sq_big.data[a * block_size_big_mul * (block_size * block_size_big_mul) + b * block_size_big_mul]; | |
| auto o = sq_small.data[a * block_size + b] - q * t - p; | |
| dist += o * o; | |
| } | |
| } | |
| if (dist < best_dist) { | |
| size_t j_transform = j % s_b; | |
| best_dist = dist; | |
| best_idx = j_transform; | |
| best_t = static_cast<transformation>(j / s_b); | |
| best_q = q; | |
| best_p = p; | |
| } | |
| } | |
| //printf("best_idx = %d\n", best_idx); | |
| res[i] = mapping(best_idx, best_t, best_q, best_p, best_dist); | |
| } | |
| return res; | |
| } | |
| auto fill_with_mappings( | |
| const std::vector<square>& par_16x16, | |
| const std::vector<square>& par_8x8, | |
| const std::vector<square>& par_4x4, | |
| const std::vector<square>& parbig_32x32, | |
| const std::vector<square>& parbig_16x16, | |
| const std::vector<square>& parbig_8x8, | |
| quadtree_node& node) | |
| { | |
| assert(node.my_block_size <= 16); | |
| assert(node.my_block_size >= 4); | |
| node.leftup_is_quadtree = false; | |
| node.rightup_is_quadtree = false; | |
| node.leftdown_is_quadtree = false; | |
| node.rightdown_is_quadtree = false; | |
| auto& par_small = (node.my_block_size == 16) ? par_16x16 : ((node.my_block_size == 8) ? par_8x8 : par_4x4); | |
| auto& par_big = (node.my_block_size == 16) ? parbig_32x32 : ((node.my_block_size == 8) ? parbig_16x16 : parbig_8x8); | |
| auto par_small_width = static_cast<size_t>(sqrt(par_small.size())); | |
| std::vector<square> par{ | |
| par_small[node.y * 2 * par_small_width + node.x * 2], | |
| par_small[node.y * 2 * par_small_width + node.x * 2 + 1], | |
| par_small[node.y * 2 * par_small_width + par_small_width + node.x * 2], | |
| par_small[node.y * 2 * par_small_width + par_small_width + node.x * 2 + 1] | |
| }; | |
| auto m = make_mapping(par, par_big, node.my_block_size, 2); | |
| node.lu.map = new mapping(std::move(m[0])); | |
| node.ru.map = new mapping(std::move(m[1])); | |
| node.ld.map = new mapping(std::move(m[2])); | |
| node.rd.map = new mapping(std::move(m[3])); | |
| } | |
| void fill( | |
| const std::vector<square>& par_16x16, | |
| const std::vector<square>& par_8x8, | |
| const std::vector<square>& par_4x4, | |
| const std::vector<square>& parbig_32x32, | |
| const std::vector<square>& parbig_16x16, | |
| const std::vector<square>& parbig_8x8, | |
| quadtree_node& node, | |
| std::vector<std::pair<mapping*, mapping_parent>>& bottom_level_nodes) | |
| { | |
| assert(node.my_block_size > 16); | |
| node.leftup_is_quadtree = true; | |
| node.rightup_is_quadtree = true; | |
| node.leftdown_is_quadtree = true; | |
| node.rightdown_is_quadtree = true; | |
| node.lu.node = new quadtree_node(node.my_block_size / 2); | |
| node.ru.node = new quadtree_node(node.my_block_size / 2); | |
| node.ld.node = new quadtree_node(node.my_block_size / 2); | |
| node.rd.node = new quadtree_node(node.my_block_size / 2); | |
| node.lu.node->x = node.x * 2; | |
| node.lu.node->y = node.y * 2; | |
| node.ru.node->x = node.x * 2 + 1; | |
| node.ru.node->y = node.y * 2; | |
| node.ld.node->x = node.x * 2; | |
| node.ld.node->y = node.y * 2 + 1; | |
| node.rd.node->x = node.x * 2 + 1; | |
| node.rd.node->y = node.y * 2 + 1; | |
| if (node.my_block_size > 32) { | |
| fill( | |
| par_16x16, | |
| par_8x8, | |
| par_4x4, | |
| parbig_32x32, | |
| parbig_16x16, | |
| parbig_8x8, | |
| *node.lu.node, | |
| bottom_level_nodes); | |
| fill( | |
| par_16x16, | |
| par_8x8, | |
| par_4x4, | |
| parbig_32x32, | |
| parbig_16x16, | |
| parbig_8x8, | |
| *node.ru.node, | |
| bottom_level_nodes); | |
| fill( | |
| par_16x16, | |
| par_8x8, | |
| par_4x4, | |
| parbig_32x32, | |
| parbig_16x16, | |
| parbig_8x8, | |
| *node.ld.node, | |
| bottom_level_nodes); | |
| fill( | |
| par_16x16, | |
| par_8x8, | |
| par_4x4, | |
| parbig_32x32, | |
| parbig_16x16, | |
| parbig_8x8, | |
| *node.rd.node, | |
| bottom_level_nodes); | |
| } else { | |
| fill_with_mappings( | |
| par_16x16, | |
| par_8x8, | |
| par_4x4, | |
| parbig_32x32, | |
| parbig_16x16, | |
| parbig_8x8, | |
| *node.lu.node); | |
| fill_with_mappings( | |
| par_16x16, | |
| par_8x8, | |
| par_4x4, | |
| parbig_32x32, | |
| parbig_16x16, | |
| parbig_8x8, | |
| *node.ru.node); | |
| fill_with_mappings( | |
| par_16x16, | |
| par_8x8, | |
| par_4x4, | |
| parbig_32x32, | |
| parbig_16x16, | |
| parbig_8x8, | |
| *node.ld.node); | |
| fill_with_mappings( | |
| par_16x16, | |
| par_8x8, | |
| par_4x4, | |
| parbig_32x32, | |
| parbig_16x16, | |
| parbig_8x8, | |
| *node.rd.node); | |
| auto f = [&](quadtree_node* n) { | |
| bottom_level_nodes.emplace_back(n->lu.map, mapping_parent(n, mapping_parent::LU)); | |
| bottom_level_nodes.emplace_back(n->ru.map, mapping_parent(n, mapping_parent::RU)); | |
| bottom_level_nodes.emplace_back(n->ld.map, mapping_parent(n, mapping_parent::LD)); | |
| bottom_level_nodes.emplace_back(n->rd.map, mapping_parent(n, mapping_parent::RD)); | |
| }; | |
| f(node.lu.node); | |
| f(node.ru.node); | |
| f(node.ld.node); | |
| f(node.rd.node); | |
| } | |
| } | |
| auto make_quadtree( | |
| const std::vector<square>& par_16x16, | |
| const std::vector<square>& par_8x8, | |
| const std::vector<square>& par_4x4, | |
| const std::vector<square>& parbig_32x32, | |
| const std::vector<square>& parbig_16x16, | |
| const std::vector<square>& parbig_8x8, | |
| size_t w, | |
| int quality) | |
| { | |
| quadtree_node top(w / 2); | |
| top.x = 0; | |
| top.y = 0; | |
| std::vector<std::pair<mapping*, mapping_parent>> bottom_nodes; | |
| fill(par_16x16, par_8x8, par_4x4, parbig_32x32, parbig_16x16, parbig_8x8, top, bottom_nodes); | |
| std::sort(bottom_nodes.begin(), bottom_nodes.end(), [](auto a, auto b) { | |
| return a.first->dist > b.first->dist; | |
| }); | |
| std::vector<std::pair<mapping*, mapping_parent>> bottom_nodes_2; | |
| for (size_t i = 0; i < (bottom_nodes.size() * quality) / 100; ++i) { | |
| auto p = bottom_nodes[i]; | |
| delete p.first; | |
| auto mappar = p.second; | |
| auto node = mappar.parent; | |
| switch (mappar.position) { | |
| case mapping_parent::LU: | |
| node->leftup_is_quadtree = true; | |
| node->lu.node = new quadtree_node(node->my_block_size / 2); | |
| node->lu.node->x = node->x * 2; | |
| node->lu.node->y = node->y * 2; | |
| fill_with_mappings( | |
| par_16x16, | |
| par_8x8, | |
| par_4x4, | |
| parbig_32x32, | |
| parbig_16x16, | |
| parbig_8x8, | |
| *node->lu.node); | |
| bottom_nodes_2.emplace_back(node->lu.node->lu.map, mapping_parent(node->lu.node, mapping_parent::LU)); | |
| bottom_nodes_2.emplace_back(node->lu.node->ru.map, mapping_parent(node->lu.node, mapping_parent::RU)); | |
| bottom_nodes_2.emplace_back(node->lu.node->ld.map, mapping_parent(node->lu.node, mapping_parent::LD)); | |
| bottom_nodes_2.emplace_back(node->lu.node->rd.map, mapping_parent(node->lu.node, mapping_parent::RD)); | |
| break; | |
| case mapping_parent::RU: | |
| node->rightup_is_quadtree = true; | |
| node->ru.node = new quadtree_node(node->my_block_size / 2); | |
| node->ru.node->x = node->x * 2 + 1; | |
| node->ru.node->y = node->y * 2; | |
| fill_with_mappings( | |
| par_16x16, | |
| par_8x8, | |
| par_4x4, | |
| parbig_32x32, | |
| parbig_16x16, | |
| parbig_8x8, | |
| *node->ru.node); | |
| bottom_nodes_2.emplace_back(node->ru.node->lu.map, mapping_parent(node->ru.node, mapping_parent::LU)); | |
| bottom_nodes_2.emplace_back(node->ru.node->ru.map, mapping_parent(node->ru.node, mapping_parent::RU)); | |
| bottom_nodes_2.emplace_back(node->ru.node->ld.map, mapping_parent(node->ru.node, mapping_parent::LD)); | |
| bottom_nodes_2.emplace_back(node->ru.node->rd.map, mapping_parent(node->ru.node, mapping_parent::RD)); | |
| break; | |
| case mapping_parent::LD: | |
| node->leftdown_is_quadtree = true; | |
| node->ld.node = new quadtree_node(node->my_block_size / 2); | |
| node->ld.node->x = node->x * 2; | |
| node->ld.node->y = node->y * 2 + 1; | |
| fill_with_mappings( | |
| par_16x16, | |
| par_8x8, | |
| par_4x4, | |
| parbig_32x32, | |
| parbig_16x16, | |
| parbig_8x8, | |
| *node->ld.node); | |
| bottom_nodes_2.emplace_back(node->ld.node->lu.map, mapping_parent(node->ld.node, mapping_parent::LU)); | |
| bottom_nodes_2.emplace_back(node->ld.node->ru.map, mapping_parent(node->ld.node, mapping_parent::RU)); | |
| bottom_nodes_2.emplace_back(node->ld.node->ld.map, mapping_parent(node->ld.node, mapping_parent::LD)); | |
| bottom_nodes_2.emplace_back(node->ld.node->rd.map, mapping_parent(node->ld.node, mapping_parent::RD)); | |
| break; | |
| case mapping_parent::RD: | |
| node->rightdown_is_quadtree = true; | |
| node->rd.node = new quadtree_node(node->my_block_size / 2); | |
| node->rd.node->x = node->x * 2 + 1; | |
| node->rd.node->y = node->y * 2 + 1; | |
| fill_with_mappings( | |
| par_16x16, | |
| par_8x8, | |
| par_4x4, | |
| parbig_32x32, | |
| parbig_16x16, | |
| parbig_8x8, | |
| *node->rd.node); | |
| bottom_nodes_2.emplace_back(node->rd.node->lu.map, mapping_parent(node->rd.node, mapping_parent::LU)); | |
| bottom_nodes_2.emplace_back(node->rd.node->ru.map, mapping_parent(node->rd.node, mapping_parent::RU)); | |
| bottom_nodes_2.emplace_back(node->rd.node->ld.map, mapping_parent(node->rd.node, mapping_parent::LD)); | |
| bottom_nodes_2.emplace_back(node->rd.node->rd.map, mapping_parent(node->rd.node, mapping_parent::RD)); | |
| break; | |
| } | |
| } | |
| std::sort(bottom_nodes_2.begin(), bottom_nodes_2.end(), [](auto a, auto b) { | |
| return a.first->dist > b.first->dist; | |
| }); | |
| for (size_t i = 0; i < (bottom_nodes_2.size() * quality) / 100; ++i) { | |
| auto p = bottom_nodes_2[i]; | |
| delete p.first; | |
| auto mappar = p.second; | |
| auto node = mappar.parent; | |
| switch (mappar.position) { | |
| case mapping_parent::LU: | |
| node->leftup_is_quadtree = true; | |
| node->lu.node = new quadtree_node(node->my_block_size / 2); | |
| node->lu.node->x = node->x * 2; | |
| node->lu.node->y = node->y * 2; | |
| fill_with_mappings( | |
| par_16x16, | |
| par_8x8, | |
| par_4x4, | |
| parbig_32x32, | |
| parbig_16x16, | |
| parbig_8x8, | |
| *node->lu.node); | |
| break; | |
| case mapping_parent::RU: | |
| node->rightup_is_quadtree = true; | |
| node->ru.node = new quadtree_node(node->my_block_size / 2); | |
| node->ru.node->x = node->x * 2 + 1; | |
| node->ru.node->y = node->y * 2; | |
| fill_with_mappings( | |
| par_16x16, | |
| par_8x8, | |
| par_4x4, | |
| parbig_32x32, | |
| parbig_16x16, | |
| parbig_8x8, | |
| *node->ru.node); | |
| break; | |
| case mapping_parent::LD: | |
| node->leftdown_is_quadtree = true; | |
| node->ld.node = new quadtree_node(node->my_block_size / 2); | |
| node->ld.node->x = node->x * 2; | |
| node->ld.node->y = node->y * 2 + 1; | |
| fill_with_mappings( | |
| par_16x16, | |
| par_8x8, | |
| par_4x4, | |
| parbig_32x32, | |
| parbig_16x16, | |
| parbig_8x8, | |
| *node->ld.node); | |
| break; | |
| case mapping_parent::RD: | |
| node->rightdown_is_quadtree = true; | |
| node->rd.node = new quadtree_node(node->my_block_size / 2); | |
| node->rd.node->x = node->x * 2 + 1; | |
| node->rd.node->y = node->y * 2 + 1; | |
| fill_with_mappings( | |
| par_16x16, | |
| par_8x8, | |
| par_4x4, | |
| parbig_32x32, | |
| parbig_16x16, | |
| parbig_8x8, | |
| *node->rd.node); | |
| break; | |
| } | |
| } | |
| return top; | |
| } | |
| void del_quadtree(quadtree_node& top) | |
| { | |
| if (top.leftup_is_quadtree) { | |
| del_quadtree(*top.lu.node); | |
| delete top.lu.node; | |
| } else { | |
| delete top.lu.map; | |
| } | |
| if (top.rightup_is_quadtree) { | |
| del_quadtree(*top.ru.node); | |
| delete top.ru.node; | |
| } else { | |
| delete top.ru.map; | |
| } | |
| if (top.leftdown_is_quadtree) { | |
| del_quadtree(*top.ld.node); | |
| delete top.ld.node; | |
| } else { | |
| delete top.ld.map; | |
| } | |
| if (top.rightdown_is_quadtree) { | |
| del_quadtree(*top.rd.node); | |
| delete top.rd.node; | |
| } else { | |
| delete top.rd.map; | |
| } | |
| } | |
| void iterate(std::vector<double>& y, const std::vector<mapping>& map, size_t w, size_t block_size, size_t block_size_big_mul, double minval, double maxval) | |
| { | |
| auto y_copy = y; | |
| auto par = partition(y_copy, w, block_size * block_size_big_mul); | |
| size_t s_m = w / block_size; | |
| #pragma omp parallel for | |
| for (int j = 0; j < static_cast<int>(s_m); ++j) { // omp needs a signed int. | |
| for (size_t i = 0; i < s_m; ++i) { | |
| auto m = map[j * s_m + i]; | |
| //printf("Block %d, %d matched from block %d, %d with transformation = %d\n", i, j, m.idx_big % (s_m / 2), m.idx_big / (s_m / 2), m.t); | |
| auto sq = par[m.idx_big].transformed(m.t); | |
| //printf("x_c = %d;\ty_c = %d;\tt = %d\n", m.x_big, m.y_big, static_cast<int>(m.t)); | |
| for (size_t y_ = 0; y_ < block_size; ++y_) { | |
| for (size_t x_ = 0; x_ < block_size; ++x_) { | |
| auto temp = m.q * sq.data[y_ * block_size_big_mul * (block_size * block_size_big_mul) + x_ * block_size_big_mul] + m.p; | |
| temp = (temp > maxval) ? maxval : temp; | |
| temp = (temp < minval) ? minval : temp; | |
| y[(j * block_size + y_) * w + (i * block_size + x_)] = temp; | |
| } | |
| } | |
| } | |
| } | |
| } | |
| auto iterate(const std::vector<double>& y, const std::vector<quadtree_node*> nodes, size_t w, size_t w_out, double minval, double maxval) | |
| { | |
| std::vector<double> output(w_out * w_out, 0.0); | |
| auto par_32x32 = partition(y, w, 32); | |
| auto par_16x16 = partition(y, w, 16); | |
| auto par_8x8 = partition(y, w, 8); | |
| #pragma omp parallel for | |
| for (int i = 0; i < static_cast<int>(nodes.size()); ++i) { // omp needs a signed int. | |
| auto node = nodes[i]; | |
| auto& par = (node->my_block_size == 16) ? par_32x32 : ((node->my_block_size == 8) ? par_16x16 : par_8x8); | |
| if (!node->leftup_is_quadtree) { | |
| auto m = node->lu.map; | |
| auto sq = par[m->idx_big].transformed(m->t); | |
| for (size_t y_ = 0; y_ < node->my_block_size; ++y_) { | |
| for (size_t x_ = 0; x_ < node->my_block_size; ++x_) { | |
| auto temp = m->q * sq.data[y_ * 2 * (node->my_block_size * 2) + x_ * 2] + m->p; | |
| temp = (temp > maxval) ? maxval : temp; | |
| temp = (temp < minval) ? minval : temp; | |
| output[(node->y * 2 * node->my_block_size + y_) * w_out + (node->x * 2 * node->my_block_size + x_)] = temp; | |
| } | |
| } | |
| } | |
| if (!node->rightup_is_quadtree) { | |
| auto m = node->ru.map; | |
| auto sq = par[m->idx_big].transformed(m->t); | |
| for (size_t y_ = 0; y_ < node->my_block_size; ++y_) { | |
| for (size_t x_ = 0; x_ < node->my_block_size; ++x_) { | |
| auto temp = m->q * sq.data[y_ * 2 * (node->my_block_size * 2) + x_ * 2] + m->p; | |
| temp = (temp > maxval) ? maxval : temp; | |
| temp = (temp < minval) ? minval : temp; | |
| output[(node->y * 2 * node->my_block_size + y_) * w_out + (node->x * 2 * node->my_block_size + node->my_block_size + x_)] = temp; | |
| } | |
| } | |
| } | |
| if (!node->leftdown_is_quadtree) { | |
| auto m = node->ld.map; | |
| auto sq = par[m->idx_big].transformed(m->t); | |
| for (size_t y_ = 0; y_ < node->my_block_size; ++y_) { | |
| for (size_t x_ = 0; x_ < node->my_block_size; ++x_) { | |
| auto temp = m->q * sq.data[y_ * 2 * (node->my_block_size * 2) + x_ * 2] + m->p; | |
| temp = (temp > maxval) ? maxval : temp; | |
| temp = (temp < minval) ? minval : temp; | |
| output[(node->y * 2 * node->my_block_size + node->my_block_size + y_) * w_out + (node->x * 2 * node->my_block_size + x_)] = temp; | |
| } | |
| } | |
| } | |
| if (!node->rightdown_is_quadtree) { | |
| auto m = node->rd.map; | |
| auto sq = par[m->idx_big].transformed(m->t); | |
| for (size_t y_ = 0; y_ < node->my_block_size; ++y_) { | |
| for (size_t x_ = 0; x_ < node->my_block_size; ++x_) { | |
| auto temp = m->q * sq.data[y_ * 2 * (node->my_block_size * 2) + x_ * 2] + m->p; | |
| temp = (temp > maxval) ? maxval : temp; | |
| temp = (temp < minval) ? minval : temp; | |
| output[(node->y * 2 * node->my_block_size + node->my_block_size + y_) * w_out + (node->x * 2 * node->my_block_size + node->my_block_size + x_)] = temp; | |
| } | |
| } | |
| } | |
| } | |
| return output; | |
| } | |
| void write_mapping(BitVec& data, const mapping& m, size_t w, size_t block_size) | |
| { | |
| size_t max_idx = (w / (block_size / 2)) * (w / (block_size / 2)); | |
| size_t bits_for_idx = std::lround(std::log2(static_cast<double>(max_idx))); | |
| for (size_t i = 0; i < bits_for_idx; ++i) | |
| data.add_bit((m.idx_big & (1 << i)) > 0); | |
| data.add_bit((m.t & 1) > 0); | |
| data.add_bit((m.t & 2) > 0); | |
| data.add_bit((m.t & 4) > 0); | |
| // 1 byte for q. | |
| auto temp = (m.q - Q_MIN) * (255 / (Q_MAX - Q_MIN)); | |
| if (temp > 255.0) | |
| temp = 255.0; | |
| else if (temp < 0.0) | |
| temp = 0.0; | |
| auto q = static_cast<uint8_t>(temp); | |
| data.add_byte(q); | |
| // 1 byte for p. | |
| double p; | |
| if (m.p >= 0) | |
| p = sqrt(m.p); | |
| else | |
| p = -sqrt(-m.p); | |
| temp = (p + 1.0) * (255 / 2.0); | |
| if (temp > 255.0) | |
| temp = 255.0; | |
| else if (temp < 0.0) | |
| temp = 0.0; | |
| data.add_byte(static_cast<uint8_t>(temp)); | |
| } | |
| void read_mapping(BitArr& data, mapping& m, size_t w, size_t block_size) | |
| { | |
| size_t max_idx = (w / (block_size / 2)) * (w / (block_size / 2)); | |
| size_t bits_for_idx = std::lround(std::log2(static_cast<double>(max_idx))); | |
| m.idx_big = 0; | |
| for (size_t i = 0; i < bits_for_idx; ++i) | |
| m.idx_big |= data.get_bit() << i; | |
| int t = data.get_bit(); | |
| t |= data.get_bit() << 1; | |
| t |= data.get_bit() << 2; | |
| m.t = static_cast<transformation>(t); | |
| uint8_t q = data.get_byte(); | |
| m.q = q * ((Q_MAX - Q_MIN) / 255) + Q_MIN; | |
| uint8_t p = data.get_byte(); | |
| auto p_ = p * (2.0 / 255) - 1.0; | |
| if (p_ >= 0) | |
| m.p = p_ * p_; | |
| else | |
| m.p = -(p_ * p_); | |
| //printf("%.8f\n", m.p); | |
| } | |
| void write_quadtree(BitVec& data, quadtree_node& top, size_t w) | |
| { | |
| data.add_bit(top.leftup_is_quadtree); | |
| data.add_bit(top.rightup_is_quadtree); | |
| data.add_bit(top.leftdown_is_quadtree); | |
| data.add_bit(top.rightdown_is_quadtree); | |
| if (top.leftup_is_quadtree) { | |
| write_quadtree(data, *top.lu.node, w); | |
| } else { | |
| write_mapping(data, *top.lu.map, w, top.my_block_size); | |
| } | |
| if (top.rightup_is_quadtree) { | |
| write_quadtree(data, *top.ru.node, w); | |
| } else { | |
| write_mapping(data, *top.ru.map, w, top.my_block_size); | |
| } | |
| if (top.leftdown_is_quadtree) { | |
| write_quadtree(data, *top.ld.node, w); | |
| } else { | |
| write_mapping(data, *top.ld.map, w, top.my_block_size); | |
| } | |
| if (top.rightdown_is_quadtree) { | |
| write_quadtree(data, *top.rd.node, w); | |
| } else { | |
| write_mapping(data, *top.rd.map, w, top.my_block_size); | |
| } | |
| } | |
| void read_quadtree(BitArr& data, size_t w, size_t cur_block_size, int cur_x, int cur_y, quadtree_node*& out, std::vector<quadtree_node*>& bottom_level_nodes) | |
| { | |
| out = new quadtree_node(cur_block_size); | |
| out->x = cur_x; | |
| out->y = cur_y; | |
| out->leftup_is_quadtree = data.get_bit(); | |
| out->rightup_is_quadtree = data.get_bit(); | |
| out->leftdown_is_quadtree = data.get_bit(); | |
| out->rightdown_is_quadtree = data.get_bit(); | |
| if (!out->leftup_is_quadtree | |
| || !out->rightup_is_quadtree | |
| || !out->leftdown_is_quadtree | |
| || !out->rightdown_is_quadtree) { | |
| bottom_level_nodes.push_back(out); | |
| } | |
| if (out->leftup_is_quadtree) { | |
| read_quadtree(data, w, cur_block_size / 2, cur_x * 2, cur_y * 2, out->lu.node, bottom_level_nodes); | |
| } else { | |
| mapping m(0, NONE, 0.0, 0.0); | |
| read_mapping(data, m, w, cur_block_size); | |
| out->lu.map = new mapping(std::move(m)); | |
| } | |
| if (out->rightup_is_quadtree) { | |
| read_quadtree(data, w, cur_block_size / 2, cur_x * 2 + 1, cur_y * 2, out->ru.node, bottom_level_nodes); | |
| } else { | |
| mapping m(0, NONE, 0.0, 0.0); | |
| read_mapping(data, m, w, cur_block_size); | |
| out->ru.map = new mapping(std::move(m)); | |
| } | |
| if (out->leftdown_is_quadtree) { | |
| read_quadtree(data, w, cur_block_size / 2, cur_x * 2, cur_y * 2 + 1, out->ld.node, bottom_level_nodes); | |
| } else { | |
| mapping m(0, NONE, 0.0, 0.0); | |
| read_mapping(data, m, w, cur_block_size); | |
| out->ld.map = new mapping(std::move(m)); | |
| } | |
| if (out->rightdown_is_quadtree) { | |
| read_quadtree(data, w, cur_block_size / 2, cur_x * 2 + 1, cur_y * 2 + 1, out->rd.node, bottom_level_nodes); | |
| } else { | |
| mapping m(0, NONE, 0.0, 0.0); | |
| read_mapping(data, m, w, cur_block_size); | |
| out->rd.map = new mapping(std::move(m)); | |
| } | |
| } | |
| /* | |
| This function should compress image data. | |
| Arguments: | |
| image - pointer to raw rgb image data. To access pixel at point x,y use image[y*w + x]. Pixel format 0x00RRGGBB | |
| compressed_file_name - open this file to write your output | |
| w - width of input image | |
| h - height of input image | |
| block_size - size of grid block | |
| quality - compresion quality. Range [1,100] | |
| */ | |
| void Compress(unsigned long* image, const char* compressed_file_name, size_t w, size_t block_size, int quality) | |
| { | |
| std::vector<double> y, u, v; | |
| bool grayscale = colorspaces::yuv_from_rgb_discr(w, w, image, y, u, v); | |
| if (quality == 100) { | |
| quality = 0; | |
| block_size = 4; | |
| } | |
| std::vector<unsigned char> output_vec; | |
| BitVec bits(output_vec); | |
| bits.add_bit(w == 512); | |
| bits.add_bit(quality != 0); | |
| bits.add_bit(grayscale); | |
| if (quality == 0) { | |
| bits.add_bit(block_size == 4 || block_size == 16); | |
| bits.add_bit(block_size == 16); | |
| auto par_small = partition(y, w, block_size); | |
| auto par_large = partition(y, w, block_size * 2); | |
| add_transforms(par_large); | |
| auto m = make_mapping(par_small, par_large, block_size, 2); | |
| size_t s = m.size(); | |
| for (size_t i = 0; i < s; ++i) { | |
| write_mapping(bits, m[i], w, block_size); | |
| } | |
| if (!grayscale) { | |
| par_small = partition(u, w / 2, block_size); | |
| m = make_mapping(par_small, par_large, block_size, 2); | |
| s = s / 4; | |
| for (size_t i = 0; i < s; ++i) { | |
| write_mapping(bits, m[i], w, block_size); | |
| } | |
| par_small = partition(v, w / 2, block_size); | |
| m = make_mapping(par_small, par_large, block_size, 2); | |
| for (size_t i = 0; i < s; ++i) { | |
| write_mapping(bits, m[i], w, block_size); | |
| } | |
| } | |
| } else { | |
| auto par_16x16 = partition(y, w, 16); | |
| auto par_8x8 = partition(y, w, 8); | |
| auto par_4x4 = partition(y, w, 4); | |
| auto parbig_32x32 = partition(y, w, 32); | |
| auto parbig_16x16 = par_16x16; | |
| auto parbig_8x8 = par_8x8; | |
| add_transforms(parbig_32x32); | |
| add_transforms(parbig_16x16); | |
| add_transforms(parbig_8x8); | |
| auto tree = make_quadtree(par_16x16, par_8x8, par_4x4, parbig_32x32, parbig_16x16, parbig_8x8, w, quality); | |
| write_quadtree(bits, tree, w); | |
| del_quadtree(tree); | |
| if (!grayscale) { | |
| par_16x16 = partition(u, w / 2, 16); | |
| par_8x8 = partition(u, w / 2, 8); | |
| par_4x4 = partition(u, w / 2, 4); | |
| tree = make_quadtree(par_16x16, par_8x8, par_4x4, parbig_32x32, parbig_16x16, parbig_8x8, w / 2, quality); | |
| write_quadtree(bits, tree, w); | |
| del_quadtree(tree); | |
| par_16x16 = partition(v, w / 2, 16); | |
| par_8x8 = partition(v, w / 2, 8); | |
| par_4x4 = partition(v, w / 2, 4); | |
| tree = make_quadtree(par_16x16, par_8x8, par_4x4, parbig_32x32, parbig_16x16, parbig_8x8, w / 2, quality); | |
| write_quadtree(bits, tree, w); | |
| del_quadtree(tree); | |
| } | |
| } | |
| bits.align(); | |
| FILE *fp = fopen(compressed_file_name, "wb"); | |
| printf("Starting arithmetic compression\n"); | |
| compress::compress(output_vec.data(), output_vec.size(), fp, false, false); | |
| fclose(fp); | |
| } | |
| /* | |
| This function should decompress image data | |
| Arguments: | |
| compressed_file_name - open this file to read your compressed data from it | |
| w - set this variable to the actual width of output image | |
| h - set this variable to the actual height of output image | |
| Return value: | |
| Pointer to raw output image data (Pixel format 0x00RRGGBB) | |
| */ | |
| unsigned long* Decompress(const char* compressed_file_name, size_t& w) | |
| { | |
| FILE *fp = fopen(compressed_file_name, "rb"); | |
| std::vector<uint8_t> in; | |
| compress::decompress(fp, in); | |
| printf("Arithmetic decompression complete\n"); | |
| fclose(fp); | |
| BitArr bits(in.data()); | |
| w = bits.get_bit() ? 512 : 256; | |
| bool quadtree = bits.get_bit(); | |
| bool grayscale = bits.get_bit(); | |
| unsigned long *image = new unsigned long[w * w]; | |
| std::vector<double> y(w * w, 0.5); | |
| if (!quadtree) { | |
| size_t s; | |
| size_t block_size; | |
| if (bits.get_bit()) { | |
| block_size = bits.get_bit() ? 16 : 4; | |
| } else { | |
| bits.get_bit(); | |
| block_size = 8; | |
| } | |
| s = (w / block_size) * (w / block_size); | |
| std::vector<mapping> m(s, mapping(0, NONE, 0.0, 0.0)); | |
| for (size_t i = 0; i < s; ++i) { | |
| read_mapping(bits, m[i], w, block_size); | |
| } | |
| for (int i = 0; i < 100; ++i) { | |
| iterate(y, m, w, block_size, 2, 0.0, 1.0); | |
| } | |
| if (grayscale) { | |
| colorspaces::rgb_from_y(w, w, y, image); | |
| } else { | |
| s = s / 4; | |
| std::vector<mapping> m_u(s, mapping(0, NONE, 0.0, 0.0)); | |
| for (size_t i = 0; i < s; ++i) { | |
| read_mapping(bits, m_u[i], w, block_size); | |
| } | |
| std::vector<mapping> m_v(s, mapping(0, NONE, 0.0, 0.0)); | |
| for (size_t i = 0; i < s; ++i) { | |
| read_mapping(bits, m_v[i], w, block_size); | |
| } | |
| //auto u = y, v = y; | |
| //iterate(u, m_u, w, block_size, 2, -0.436, 0.436); | |
| //iterate(v, m_v, w, block_size, 2, -0.615, 0.615); | |
| std::vector<double> u(w * w / 4); | |
| std::vector<double> v(w * w / 4); | |
| auto y_copy = y; | |
| auto par = partition(y_copy, w, block_size * 2); | |
| size_t s_m = (w / 2) / block_size; | |
| #pragma omp parallel for | |
| for (int j = 0; j < static_cast<int>(s_m); ++j) { // omp needs a signed int. | |
| for (size_t i = 0; i < s_m; ++i) { | |
| auto m = m_u[j * s_m + i]; | |
| auto m2 = m_v[j * s_m + i]; | |
| auto sq = par[m.idx_big].transformed(m.t); | |
| auto sq2 = par[m2.idx_big].transformed(m2.t); | |
| //printf("x_c = %d;\ty_c = %d;\tt = %d\n", m.x_big, m.y_big, static_cast<int>(m.t)); | |
| for (size_t y_ = 0; y_ < block_size; ++y_) { | |
| for (size_t x_ = 0; x_ < block_size; ++x_) { | |
| auto temp = m.q * sq.data[y_ * 2 * (block_size * 2) + x_ * 2] + m.p; | |
| auto temp2 = m2.q * sq2.data[y_ * 2 * (block_size * 2) + x_ * 2] + m2.p; | |
| temp = (temp > 0.436) ? 0.436 : temp; | |
| temp = (temp < -0.436) ? -0.436 : temp; | |
| temp2 = (temp2 > 0.615) ? 0.615 : temp2; | |
| temp2 = (temp2 < -0.615) ? -0.615 : temp2; | |
| u[(j * block_size + y_) * (w / 2) + (i * block_size + x_)] = temp; | |
| v[(j * block_size + y_) * (w / 2) + (i * block_size + x_)] = temp2; | |
| } | |
| } | |
| } | |
| } | |
| colorspaces::rgb_from_yuv_discr(w, w, y, u, v, image); | |
| } | |
| } else { | |
| quadtree_node *tree; | |
| std::vector<quadtree_node*> bottom_level_nodes; | |
| read_quadtree(bits, w, w / 2, 0, 0, tree, bottom_level_nodes); | |
| for (int i = 0; i < 100; ++i) { | |
| y = iterate(y, bottom_level_nodes, w, w, 0.0, 1.0); | |
| } | |
| del_quadtree(*tree); | |
| delete tree; | |
| if (grayscale) { | |
| colorspaces::rgb_from_y(w, w, y, image); | |
| } else { | |
| bottom_level_nodes.clear(); | |
| read_quadtree(bits, w, w / 4, 0, 0, tree, bottom_level_nodes); | |
| auto u = iterate(y, bottom_level_nodes, w, w / 2, -0.436, 0.436); | |
| del_quadtree(*tree); | |
| delete tree; | |
| bottom_level_nodes.clear(); | |
| read_quadtree(bits, w, w / 4, 0, 0, tree, bottom_level_nodes); | |
| auto v = iterate(y, bottom_level_nodes, w, w / 2, -0.615, 0.615); | |
| del_quadtree(*tree); | |
| delete tree; | |
| colorspaces::rgb_from_yuv_discr(w, w, y, u, v, image); | |
| } | |
| } | |
| return image; | |
| } | |
| void PSNR(unsigned long *image, unsigned long *image2, int w, int h) | |
| { | |
| std::vector<double> y1; | |
| std::vector<double> y2; | |
| colorspaces::y_from_rgb(w, h, image, y1); | |
| colorspaces::y_from_rgb(w, h, image2, y2); | |
| auto y_psnr = psnr::psnr(y1, y2); | |
| printf("%.8f\n", y_psnr); | |
| } |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment