From 249abc450abebab6c112ec3d00fdcaea66e58d27 Mon Sep 17 00:00:00 2001 From: bw Date: Feb 05 2016 17:57:33 +0000 Subject: iir_blur --- diff --git a/c++/iir_blur/Makefile b/c++/iir_blur/Makefile new file mode 100644 index 0000000..bc4d738 --- /dev/null +++ b/c++/iir_blur/Makefile @@ -0,0 +1,26 @@ +DEPLIBS = + +CXXFLAGS = -O3 -g -Wall -fmessage-length=0 `pkg-config --cflags $(DEPLIBS)` + +HEADERS = \ + surface.h + +SOURCES = \ + iir_blur.cpp \ + surface.cpp + +OBJS = \ + iir_blur.o \ + surface.o + +DEPS = $(HEADERS) $(SOURCES) +LIBS = `pkg-config --libs $(DEPLIBS)` +TARGET = iir_blur + +$(TARGET): $(OBJS) + $(CXX) -o $(TARGET) $(OBJS) $(LIBS) + +all: $(TARGET) + +clean: + rm -f $(OBJS) $(TARGET) diff --git a/c++/iir_blur/iir_blur b/c++/iir_blur/iir_blur new file mode 100755 index 0000000..ad7d08a Binary files /dev/null and b/c++/iir_blur/iir_blur differ diff --git a/c++/iir_blur/iir_blur.cpp b/c++/iir_blur/iir_blur.cpp new file mode 100644 index 0000000..be45b26 --- /dev/null +++ b/c++/iir_blur/iir_blur.cpp @@ -0,0 +1,498 @@ +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include "surface.h" + + +using namespace std; + + +#define COMPARE_HELPER(a) a < other.a ? true : other.a < a ? false : + +#define COMPARE(a) a < other.a; +#define COMPARE2(a, b) COMPARE_HELPER(a) COMPARE(b) +#define COMPARE3(a, b, c) COMPARE_HELPER(a) COMPARE2(b, c) +#define COMPARE4(a, b, c, d) COMPARE_HELPER(a) COMPARE3(b, c, d) + + +const double pi = 3.14159265358979323846; +const char logfile[] = "results/coef.log"; + +string strprintf(const char *format, ...) +{ + va_list args; + va_start(args,format); + static char buffer[1024*1024]; + vsprintf(buffer, format, args); + return buffer; +} + +template +inline bool less(const T &a, const T &b, bool def) + { return a < b ? true : b < a ? false : def; } + + +double gauss(double x, double radius) { + static const double k = 1.0/sqrt(2.0*pi); + return k/radius*exp(-0.5*x*x/(radius*radius)); +} + + +class Params { +public: + union { + struct { double k0, k1, k2; }; + double k[3]; + }; + + Params(): k0(), k1(), k2() { } + Params(double k0, double k1, double k2): k0(k0), k1(k1), k2(k2) { } + + + bool operator< (const Params &other) const + { return COMPARE3(k0, k1, k2); } +}; + + +class Checker { +private: + vector blank; + vector tmp; + +public: + int maxi; + double maxv; + bool valid; + double radius; + + Checker(double radius): maxi(), maxv(), valid(), radius(radius) { + blank.resize((int)ceil(radius*10)); + tmp.resize(blank.size()); + for(int i = 0; i < (int)blank.size(); ++i) + blank[i] = gauss(i, radius); + } + + Params convert_params(const Params ¶ms) { + double a = 1.0/params.k0; + double b = a*cos(pi*params.k1); + double c = 1.0/params.k2; + + Params p; + p.k0 = -(a*a + 2.0*c*b)/(c*a*a); + p.k1 = (c + 2.0*b)/(c*a*a); + p.k2 = -1.0/(c*a*a); + + return p; + } + + double check2(const Params ¶ms) { + if ( abs(params.k0) > 0.99999999 + || abs(params.k1) > 0.99999999 + || abs(params.k2) > 0.99999999 + || params.k0 < 1e-6 + || abs(params.k2) < 1e-6 ) return INFINITY; + return check(convert_params(params), true); + } + + double check(const Params ¶ms, bool mode3a = true) { + valid = false; + + //if (mode3a && fabs(params.k0 + (params.k1 + params.k2)/(params.k0 + params.k1 + params.k2)) >= 0.999999999) + // return INFINITY; + if (!mode3a && fabs(params.k1 + params.k2/(params.k1 + params.k2)) >= 0.999999999) + return INFINITY; + if (fabs(params.k0) < 1e-10 && fabs(params.k1) < 1e-10 && fabs(params.k2) < 1e-10) + return INFINITY; + + maxv = 1.0 - params.k0 - params.k1 - params.k2; + int offset = 3; + maxi = offset; + tmp[offset] = 1.0; + for(int i = offset; i < (int)tmp.size(); ++i) { + if (mode3a) + tmp[i] = (i == offset ? maxv : 0.0) + + params.k0*tmp[i-1] + + params.k1*tmp[i-2] + + params.k2*tmp[i-3]; + else + tmp[i] = (i == offset ? 1.0 : 0.0) + + (i-1 == offset ? params.k0 : 0.0) + + params.k1*tmp[i-1] + + params.k2*tmp[i-2]; + if (i < (int)tmp.size()/2 && maxv < tmp[i]) { maxv = tmp[i]; maxi = i; } + } + + //if (maxv < 1e-20 || isnan(maxv)) return INFINITY; + //maxi = (int)round(radius) + offset; + + //double k = blank[0]/maxv; + double delta = 0.0; + for(int i = maxi; i < (int)tmp.size(); ++i) { + double t = tmp[i]; + //double d = blank[i - maxi] - t; + //delta += d*d; + //if (t < 0) delta += 5*d*d; + delta += fabs(blank[i - maxi] - t); + if (t < 0) delta += fabs(t); + } + + if (isinf(delta) || isnan(delta)) return INFINITY; + + valid = true; + return delta; + } + + + void graph(const string &filename) { + int count = blank.size(); + int scale = 256/count + 1; + + int pad = max(16, count/8); + int gw = count*scale; + int gh = gw; + + Surface::Color axis = Surface::convert_color(0, 0, 0, 1); + Surface::Color ca = Surface::convert_color(0, 0, 1, 0.75); + Surface::Color cb = Surface::convert_color(1, 0, 0, 0.75); + + Surface s(gw + 2*pad, gh + 2*pad); + s.clear(1.0, 1.0, 1.0, 1.0); + + s.set_color(axis); + s.move_to(pad, pad); + s.line_by(0, gh); + s.line_by(gw, 0); + + s.set_color(ca); + s.move_to(pad, pad); + for(int i = 0; i < count; ++i) + s.line_to(pad + i*scale, pad + (int)round((1.0 - blank[i]/blank[0])*gh)); + + if (valid) { + s.set_color(cb); + s.move_to(pad, pad); + for(int i = 0; i < count - maxi; ++i) + s.line_to(pad + i*scale, pad + (int)round((1.0 - tmp[i + maxi]/maxv)*gh)); + } + + s.save(filename); + } +}; + + +class Finder { +public: + class Entry { + public: + double value; + Params min; + Params current; + Params max; + + Finder *finder; + Entry *parent; + int level; + int index; + int max_variants; + set variants; + + Entry(): value(INFINITY), finder(), parent(), level(), index(), max_variants() { } + + bool operator< (const Entry &other) const + { return COMPARE4(value, min, current, max); } + + void add_variant(const Entry &entry) { + if (isinf(entry.value) || isnan(entry.value) || max_variants == 0) return; + if ((int)variants.size() < max_variants || entry.value < variants.rbegin()->value) { + variants.insert(entry); + if (max_variants < (int)variants.size()) { + set::iterator i = variants.end(); + variants.erase(--i); + } + } + } + + void report_index() { + if (parent) { + parent->report_index(); + cout << "[" << index << "]"; + } + } + + void report(bool open) { + if (level < finder->max_report_level) { + cout << "level " << level << " "; + report_index(); + cout << (open ? " begin" : " end") + << " " << current.k0 + << " " << current.k1 + << " " << current.k2 + << " " << value + << endl; + } + } + + double find() { + if (max_variants <= 0) return value; + + report(true); + + int i[3]; + for(i[0] = 0; i[0] < finder->sub_division; ++i[0]) + for(i[1] = 0; i[1] < finder->sub_division; ++i[1]) + for(i[2] = 0; i[2] < finder->sub_division; ++i[2]) { + Entry sub; + sub.finder = finder; + sub.parent = this; + sub.level = level + 1; + sub.max_variants = max_variants/4; + if (sub.max_variants < 1) sub.max_variants = 1; + for(int j = 0; j < 3; ++j) { + sub.min.k[j] = (max.k[j] - min.k[j])/finder->sub_division*(i[j] + 0.5 - 4.0) + min.k[j]; + sub.current.k[j] = (max.k[j] - min.k[j])/finder->sub_division*(i[j] + 0.5 + 0.0) + min.k[j]; + sub.max.k[j] = (max.k[j] - min.k[j])/finder->sub_division*(i[j] + 0.5 + 4.0) + min.k[j]; + } + sub.value = finder->checker.check2(sub.current); + add_variant(sub); + } + + if (level < finder->max_level) { + vector v(variants.begin(), variants.end()); + variants.clear(); + for(vector::iterator i = v.begin(); i != v.end(); ++i) { + i->index = i - v.begin(); + i->find(); + variants.insert(*i); + } + } + + if (!variants.empty() && variants.begin()->value < value) { + value = variants.begin()->value; + current = variants.begin()->current; + } + + report(false); + return value; + } + }; + + Checker checker; + + int sub_division; + int max_level; + int max_report_level; + + Entry root; + + Finder( + double radius, + Params min, + Params max, + int sub_division, + int max_variants, + int max_level, + int max_report_level + ): + checker(radius), + sub_division(sub_division), + max_level(max_level), + max_report_level(max_report_level) + { + root.finder = this; + root.min = min; + root.max = max; + root.max_variants = max_variants; + } + + double find() { + for(int j = 0; j < 3; ++j) + root.current.k[j] = 0.5*(root.max.k[j] - root.min.k[j]) + root.min.k[j]; + root.value = checker.check2(root.current); + return root.find(); + } +}; + +double find_near(Checker &checker, double step, Params ¶ms) { + double best_value = checker.check2(params); + Params best_params = params; + while(true) { + bool found = false; + for(int i = 0; i < 3; ++i) { + for(int j = -1; j <= 1; j += 2) { + Params p = params; + p.k[i] += j*step; + double v = checker.check2(p); + if (v < best_value) { best_value = v; best_params = p; } + } + } + if (found) params = best_params; else break; + } + return best_value; +} + +void log_begin() { + ofstream f(logfile, ios_base::app); + f << endl << "double coefs[][7] = {" << endl; +} + +void log_params(double radius, int maxi, double maxv, const Params ¶ms, double value) { + ofstream f(logfile, ios_base::app); + f << " { " << radius + << ", " << maxi + << setprecision(13) + << ", " << maxv + << ", " << params.k0 + << ", " << params.k1 + << ", " << params.k2 + << ", " << value + << " }," << endl; + cout << endl; +} + +void log_params(const Checker &checker, const Params ¶ms, double value) { + log_params(checker.radius, checker.maxi, checker.maxv, params, value); +} + +void log_params(const Finder &finder) { + log_params(finder.checker, finder.root.current, finder.root.value); +} + +void log_end() { + ofstream f(logfile, ios_base::app); + f << "};" << endl << endl; +} + +void graph(Checker &checker, const Params ¶ms) { + int ri = roundf(checker.radius); + if (fabs(checker.radius - ri) > 1e-5) return; + + bool is_power_of_two = false; + for(int i = 1; i <= 4096; i *= 2) + if (ri == i) is_power_of_two = true; + + if (!is_power_of_two && ri%16 != 0) + return; + + checker.check2(params); + checker.graph(strprintf("iir_%04d.tga", ri)); +} + +void process() { + const double initial_radius = 8.0; + //const double min_radius = 1.0; + //const double max_radius = 2048; + //const double step = 0.1; + //const double params_step = 1e-7; + + cout << "find initial params" << endl; + Finder finder( + initial_radius, + Params( 0.25, -1, -1), + Params( 1, 1, 1), + 256, + 1, + 3, + 5 + ); + + finder.find(); + finder.checker.check2(finder.root.current); + graph(finder.checker, finder.root.current); + + log_begin(); + log_params(finder); + + /* + cout << "walk back" << endl; + Params params = finder.root.current; + for(double r = finder.checker.radius; r > min_radius - 0.5*step; r -= step) { + Checker checker(r); + cout << " " << r << flush; + + double value = find_near(checker, params_step, params); + graph(checker, params); + + cout << " " << params.k0 + << " " << params.k1 + << " " << params.k2 + << " " << value + << endl; + } + + cout << "walk forward" << endl; + params = finder.root.current; + for(double r = finder.checker.radius; r < max_radius + 0.5*step; r += step) { + Checker checker(r); + cout << " " << r << flush; + + double value = find_near(checker, params_step, params); + graph(checker, params); + + cout << " " << params.k0 + << " " << params.k1 + << " " << params.k2 + << " " << value + << endl; + } + */ + + log_end(); +} + + + +int main() { + process(); + + /* + string logfile = "results/coef.log"; + + { ofstream f(logfile.c_str(), ios_base::app); f << endl << "double coefs[][7] = {" << endl; } + + Params prev(0.914, 0.129, -0.216); + for(double r = 4.0; r < 2048.1; r *= pow(2.0, 0.125)) { + cout << endl; + cout << "find coefficients for radius " << r << endl; + + double x = 0.1*log2(r)/2; + + Finder finder( + r, + Params(prev.k0 - x, prev.k1 - x, prev.k2 - x), + Params(prev.k0 + x, prev.k1 + x, prev.k2 + x), + 128, + 1, + 3, + 5 + ); + + finder.find(); + finder.checker.check(finder.root.current); + finder.checker.graph(strprintf("iir_%08.3f.tga", finder.checker.radius)); + + ofstream f(logfile.c_str(), ios_base::app); + f << " { " << finder.checker.radius + << ", " << finder.checker.maxi + << setprecision(20) + << ", " << finder.checker.maxv + << ", " << finder.root.current.k0 + << ", " << finder.root.current.k1 + << ", " << finder.root.current.k2 + << ", " << finder.root.value + << " }," << endl; + cout << endl; + + prev = finder.root.current; + } + + { ofstream f(logfile.c_str(), ios_base::app); f << "};" << endl << endl; } + */ +} diff --git a/c++/iir_blur/results/placeholder b/c++/iir_blur/results/placeholder new file mode 100644 index 0000000..e69de29 --- /dev/null +++ b/c++/iir_blur/results/placeholder diff --git a/c++/iir_blur/surface.cpp b/c++/iir_blur/surface.cpp new file mode 100644 index 0000000..54f6b0c --- /dev/null +++ b/c++/iir_blur/surface.cpp @@ -0,0 +1,127 @@ +/* + ......... 2015 Ivan Mahonin + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . +*/ + +#include +#include + +#include + +#include "surface.h" + + +using namespace std; + + +Surface::Surface(int width, int height): + buffer(new Color[width*height]), + width(width), + height(height), + color(), + x(), + y(), + blending(true) +{ + memset(buffer, 0, sizeof(Color)*width*height); +} + +Surface::~Surface() + { delete[] buffer; } + +unsigned int Surface::convert_channel(double c) { + int i = (int)round(c*255.0); + return i < 0 ? 0 : i > 255 ? 255 : (unsigned char)i; +} + +unsigned int Surface::blend_channel(unsigned int a, unsigned int dest, unsigned int src) +{ + return (dest*(255-a) + a*src) >> 8; +} + +Surface::Color Surface::convert_color(double r, double g, double b, double a) { + return convert_channel(r) + | (convert_channel(g) << 8) + | (convert_channel(b) << 16) + | (convert_channel(a) << 24); +} + +Surface::Color Surface::blend(Color dest, Color src) { + unsigned int sa = src >> 24; + if (sa >= 255) return src; else if (sa == 0) return dest; + unsigned int da = (dest >> 24) + sa; + if (da > 255) da = 255; + return blend_channel(sa, dest | 0xff, src | 0xff) + | (blend_channel(sa, (dest >> 8) | 0xff, (src >> 8) | 0xff) << 8) + | (blend_channel(sa, (dest >> 16) | 0xff, (src >> 16) | 0xff) << 16) + | (da << 24); +} + +void Surface::line(Color color, int x0, int y0, int x1, int y1, bool blending) { + if (x1 < x0) swap(x1, x0); + if (y1 < y0) swap(y1, y0); + int pdx = y1 - y0; + int pdy = x0 - x1; + int d = 0; + for(int x = x0, y = y0; x <= x1 && y <= y1;) { + if (blending) blend_point(x, y, color); else set_point(x, y, color); + if (abs(d+pdx) < abs(d+pdy)) + { d += pdx; ++x; } else { d += pdy; ++y; } + } +} + +void Surface::rect(Color color, int x0, int y0, int x1, int y1, bool blending) { + if (x1 < x0) swap(x1, x0); + if (y1 < y0) swap(y1, y0); + for(int x = x0; x < x1; ++x) + for(int y = y0; y < y1; ++y) + if (blending) blend_point(x, y, color); else set_point(x, y, color); +} + +void Surface::save(const string &filename) const +{ + // create file + ofstream f(("results/" + filename).c_str(), ofstream::out | ofstream::trunc | ofstream::binary); + + // write header + unsigned char targa_header[] = { + 0, // Length of the image ID field (0 - no ID field) + 0, // Whether a color map is included (0 - no colormap) + 2, // Compression and color types (2 - uncompressed true-color image) + 0, 0, 0, 0, 0, // Color map specification (not need for us) + 0, 0, // X-origin + 0, 0, // Y-origin + (unsigned char)(width & 0xff), // Image width + (unsigned char)(width >> 8), + (unsigned char)(height & 0xff), // Image height + (unsigned char)(height >> 8), + 32, // Bits per pixel + 0 // Image descriptor (keep zero for capability) + }; + f.write((char*)targa_header, sizeof(targa_header)); + + // write data + if (true) { + int line_size = 4*width; + const char *end = (char*)buffer; + const char *current = end + height*line_size; + while(current > end) { + current -= line_size; + f.write(current, line_size); + } + } else { + f.write((const char*)buffer, 4*width*height); + } +} diff --git a/c++/iir_blur/surface.h b/c++/iir_blur/surface.h new file mode 100644 index 0000000..541dbf5 --- /dev/null +++ b/c++/iir_blur/surface.h @@ -0,0 +1,85 @@ +/* + ......... 2016 Ivan Mahonin + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . +*/ + +#ifndef _SURFACE_H_ +#define _SURFACE_H_ + +#include + +#include + +class Surface { +public: + typedef unsigned int Color; + + Color* const buffer; + const int width; + const int height; + + Color color; + int x, y; + bool blending; + + Surface(int width, int height); + ~Surface(); + + static unsigned int convert_channel(double c); + static unsigned int blend_channel(unsigned int a, unsigned int dest, unsigned int src); + static Color convert_color(double r, double g, double b, double a = 1.0); + static Color blend(Color dest, Color src); + + bool is_inside(int x, int y) const + { return x >= 0 && x < width && y >= 0 && y < height; } + + Color& point(int x, int y) + { assert(is_inside(x, y)); return buffer[y*width + x]; } + const Color& point(int x, int y) const + { assert(is_inside(x, y)); return buffer[y*width + x]; } + + void set_color(Color color) + { this->color = color; } + void set_color(double r, double g, double b, double a = 1.0) + { set_color(convert_color(r, g, b, a)); } + void set_blending(bool blending) { this->blending = blending; } + + void line(Color color, int x0, int y0, int x1, int y1, bool blending = true); + void rect(Color color, int x0, int y0, int x1, int y1, bool blending = true); + + void move_to(int x, int y) { this->x = x; this->y = y; } + void line_to(int x, int y) { line(color, this->x, this->y, x, y, blending); move_to(x, y); }; + void rect_to(int x, int y) { rect(color, this->x, this->y, x, y, blending); move_to(x, y); }; + + void move_by(int x, int y) { move_to(this->x + x, this->y + y); } + void line_by(int x, int y) { line_to(this->x + x, this->y + y); } + void rect_by(int x, int y) { rect_to(this->x + x, this->y + y); } + + void clear(Color color = 0) { rect(color, 0, 0, width, height, false); } + void clear(double r, double g, double b, double a) { clear(convert_color(r, g, b, a)); } + + void set_point(int x, int y, unsigned int color) + { if (is_inside(x, y)) point(x, y) = color; } + void blend_point(int x, int y, unsigned int color) { + if (is_inside(x, y)) { + Color& p = point(x, y); + p = blend(p, color); + } + } + + void save(const std::string &filename) const; +}; + +#endif