summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorTom Smeding <tom.smeding@gmail.com>2020-04-28 22:10:06 +0200
committerTom Smeding <tom.smeding@gmail.com>2020-04-28 22:11:05 +0200
commitb835e0c8cc71626b10b9d0f02d9530460f82ca5e (patch)
tree417f724c4c69108bfa973c76f6b5df1b547851ea
Initial commitHEADmaster
-rw-r--r--.gitignore4
-rw-r--r--.gitmodules3
-rw-r--r--Makefile30
m---------lodepng0
-rw-r--r--main.cpp512
5 files changed, 549 insertions, 0 deletions
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..130ccff
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,4 @@
+obj/
+data.txt
+out.png
+mpmandel
diff --git a/.gitmodules b/.gitmodules
new file mode 100644
index 0000000..e1e442c
--- /dev/null
+++ b/.gitmodules
@@ -0,0 +1,3 @@
+[submodule "lodepng"]
+ path = lodepng
+ url = https://github.com/lvandeve/lodepng
diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..716021a
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,30 @@
+CXX = g++
+CXXFLAGS = -Wall -Wextra -std=c++17 -O2 -g -pthread -flto $(pkg-config --cflags mpfr)
+LDFLAGS = $(CXXFLAGS) $(pkg-config --libs mpfr)
+TARGET = mpmandel
+
+SOURCES = $(wildcard *.cpp) lodepng/lodepng.cpp
+
+OBJDIR = obj
+
+.PHONY: all clean
+
+all: $(TARGET)
+
+clean:
+ @echo "Cleaning"
+ @rm -f $(TARGET)
+ @rm -rf $(OBJDIR)
+
+
+$(OBJDIR)/%.o: %.cpp $(wildcard *.h) | $(OBJDIR)
+ @mkdir -p $(dir $@)
+ @echo "CXX $<"
+ @$(CXX) $(CXXFLAGS) -c -o $@ $<
+
+$(TARGET): $(patsubst %.cpp,$(OBJDIR)/%.o,$(SOURCES)) | $(OBJDIR)
+ @echo "LD -o $@"
+ @$(CXX) -o $@ $^ $(LDFLAGS)
+
+$(OBJDIR):
+ @mkdir -p $(OBJDIR)
diff --git a/lodepng b/lodepng
new file mode 160000
+Subproject e34ac04553e51a6982ae234d98ce6b76dd57a6a
diff --git a/main.cpp b/main.cpp
new file mode 100644
index 0000000..b3c0b64
--- /dev/null
+++ b/main.cpp
@@ -0,0 +1,512 @@
+#include <iostream>
+#include <fstream>
+#include <vector>
+#include <array>
+#include <thread>
+#include <atomic>
+#include <cstdint>
+#include <cstring>
+#include <cstdlib>
+#include <cmath>
+#include <mpfr.h>
+#include "lodepng/lodepng.h"
+
+
+#define MULTITHREAD
+
+
+constexpr size_t precision = 128;
+constexpr mpfr_rnd_t rounding = MPFR_RNDN;
+
+
+struct colorscheme_params {
+ double speed = 1.0;
+ double offset = 0.0;
+};
+
+static std::array<uint8_t, 3> colorscheme(double value, const colorscheme_params &params) {
+ if (value == -1) return {0, 0, 0};
+ else {
+ // 720 is a full circle
+ const double hue = params.offset + params.speed * value / 2;
+
+ // https://en.wikipedia.org/wiki/HSL_and_HSV#HSL_to_RGB_alternative
+ // saturation = 1.0
+ // lightness = 0.5
+ // a = saturation * min(lightness, 1 - lightness) = 0.5
+ const auto f = [hue](int n) {
+ // k = (n + hue / 30°) mod 12
+ // f(n) = lightness - a * max(-1, min(k - 3, 9 - k, 1))
+ double k = std::fmod(n + hue / 30, 12);
+ double res = 0.5 - 0.5 * std::max(-1.0, std::min(1.0, std::min(k - 3, 9 - k)));
+
+ return (uint8_t)(res * 127.99 + 127.99);
+ };
+
+ return {f(0), f(8), f(4)};
+ }
+}
+
+struct window {
+ double xmin, xmax, ymin, ymax;
+
+ window enlarge_to_ratio(double width_over_height) const {
+ window res = *this;
+
+ double current_rat = (xmax - xmin) / (ymax - ymin);
+ if (current_rat > width_over_height) {
+ double h = (xmax - xmin) / width_over_height;
+ double ymid = (ymin + ymax) / 2;
+ res.ymin = ymid - h / 2;
+ res.ymax = ymid + h / 2;
+ } else if (current_rat < width_over_height) {
+ double w = (ymax - ymin) * width_over_height;
+ double xmid = (xmin + xmax) / 2;
+ res.xmin = xmid - w / 2;
+ res.xmax = xmid + w / 2;
+ }
+
+ return res;
+ }
+
+ window recenter_x(double x) const {
+ window res = *this;
+ res.xmin = x - (xmax - xmin) / 2;
+ res.xmax = x + (xmax - xmin) / 2;
+ return res;
+ }
+
+ window recenter_y(double y) const {
+ window res = *this;
+ res.ymin = y - (ymax - ymin) / 2;
+ res.ymax = y + (ymax - ymin) / 2;
+ return res;
+ }
+
+ window set_width(double w) const {
+ window res = *this;
+ res.xmin = (xmax + xmin) / 2 - w / 2;
+ res.xmax = (xmax + xmin) / 2 + w / 2;
+ return res;
+ }
+
+ window set_height(double h) const {
+ window res = *this;
+ res.ymin = (ymax + ymin) / 2 - h / 2;
+ res.ymax = (ymax + ymin) / 2 + h / 2;
+ return res;
+ }
+};
+
+class mpfr_class {
+public:
+ mpfr_class() {
+ mpfr_init2(num, precision);
+ }
+
+ ~mpfr_class() {
+ mpfr_clear(num);
+ }
+
+ void set(const mpfr_class &a) {
+ mpfr_set(num, a.num, rounding);
+ }
+
+ void set_d(double a) {
+ mpfr_set_d(num, a, rounding);
+ }
+
+ void add(const mpfr_class &a, const mpfr_class &b) {
+ mpfr_add(num, a.num, b.num, rounding);
+ }
+
+ void sub(const mpfr_class &a, const mpfr_class &b) {
+ mpfr_sub(num, a.num, b.num, rounding);
+ }
+
+ void mul(const mpfr_class &a, const mpfr_class &b) {
+ mpfr_mul(num, a.num, b.num, rounding);
+ }
+
+ void mul_d(const mpfr_class &a, double b) {
+ mpfr_mul_d(num, a.num, b, rounding);
+ }
+
+ void fma(const mpfr_class &a, const mpfr_class &b, const mpfr_class &c) {
+ mpfr_fma(num, a.num, b.num, c.num, rounding);
+ }
+
+ void sqrt(const mpfr_class &a) {
+ mpfr_sqrt(num, a.num, rounding);
+ }
+
+ void log(const mpfr_class &a) {
+ mpfr_log(num, a.num, rounding);
+ }
+
+ void log2(const mpfr_class &a) {
+ mpfr_log2(num, a.num, rounding);
+ }
+
+ int cmp_ui(unsigned long int a) const {
+ return mpfr_cmp_ui(num, a);
+ }
+
+ bool greater(const mpfr_class &a) const {
+ return mpfr_greater_p(num, a.num);
+ }
+
+ double get_d() const {
+ return mpfr_get_d(num, rounding);
+ }
+
+private:
+ mpfr_t num;
+};
+
+template <typename T>
+class mpfr_float_shim {
+public:
+ void set(const mpfr_float_shim &a) {
+ num = a.num;
+ }
+
+ void set_d(double a) {
+ num = a;
+ }
+
+ void add(const mpfr_float_shim &a, const mpfr_float_shim &b) {
+ num = a.num + b.num;
+ }
+
+ void sub(const mpfr_float_shim &a, const mpfr_float_shim &b) {
+ num = a.num - b.num;
+ }
+
+ void mul(const mpfr_float_shim &a, const mpfr_float_shim &b) {
+ num = a.num * b.num;
+ }
+
+ void mul_d(const mpfr_float_shim &a, double b) {
+ num = a.num * b;
+ }
+
+ void fma(const mpfr_float_shim &a, const mpfr_float_shim &b, const mpfr_float_shim &c) {
+ num = a.num * b.num + c.num;
+ }
+
+ void sqrt(const mpfr_float_shim &a) {
+ num = std::sqrt(a.num);
+ }
+
+ void log(const mpfr_float_shim &a) {
+ num = std::log(a.num);
+ }
+
+ void log2(const mpfr_float_shim &a) {
+ num = std::log2(a.num);
+ }
+
+ bool greater(const mpfr_float_shim &a) const {
+ return num > a.num;
+ }
+
+ int cmp_ui(unsigned long int a) const {
+ if (num < a) return -1;
+ if (num > a) return 1;
+ return 0;
+ }
+
+ double get_d() const {
+ return num;
+ }
+
+private:
+ T num;
+};
+
+using FLT = mpfr_float_shim<double>;
+
+// Destroys the values in both escape_x and escape_y
+static double mandel_smooth_iters(size_t iters, FLT escape_x, FLT escape_y, int log_log_bound) {
+ escape_y.mul(escape_y, escape_y);
+ escape_x.fma(escape_x, escape_x, escape_y);
+ escape_x.log2(escape_x);
+ escape_x.log2(escape_x);
+ return iters + 1 + log_log_bound - escape_x.get_d();
+}
+
+struct params {
+ std::vector<float> *image;
+ size_t width, height;
+ size_t maxiter;
+ const FLT &xmin, &xmax, &ymin, &ymax;
+};
+
+static void mandel_job(const params &params, size_t from_idx, size_t to_idx) {
+ FLT x, y, a, b, a2, b2, magnitude;
+
+ // bound for smooth iteration count
+ // (source: https://www.iquilezles.org/www/articles/mset_smooth/mset_smooth.htm)
+ FLT bound_squared;
+ int log_log_bound = 3;
+ bound_squared.set_d(1 << (2 << log_log_bound)); // (2^(2^k))^2 = 2^(2^k + 2^k) = 2^(2 * 2^k)
+
+ for (size_t index = from_idx; index < to_idx; index++) {
+ const size_t xi = index % params.width;
+ const size_t yi = index / params.width;
+
+ x.sub(params.xmax, params.xmin);
+ x.mul_d(x, (double)xi / (params.width - 1));
+ x.add(x, params.xmin);
+
+ // min and max reversed here because the vertical axis is mirrored
+ y.sub(params.ymin, params.ymax);
+ y.mul_d(y, (double)yi / (params.height - 1));
+ y.add(y, params.ymax);
+
+ a.set(x);
+ b.set(y);
+ a2.mul(a, a);
+ b2.mul(b, b);
+
+ size_t iter;
+ for (iter = 0; iter < params.maxiter; iter++) {
+ b.mul(a, b);
+ b.add(b, b);
+ b.add(b, y);
+
+ a.sub(a2, b2);
+ a.add(a, x);
+
+ a2.mul(a, a);
+ b2.mul(b, b);
+
+ magnitude.add(a2, b2);
+
+ if (magnitude.greater(bound_squared)) break;
+ }
+
+ const double smooth_iters =
+ iter == params.maxiter
+ ? -1
+ : mandel_smooth_iters(iter, a, b, log_log_bound);
+
+ (*params.image)[index] = smooth_iters;
+ }
+}
+
+static std::vector<uint8_t> colorise_image(const std::vector<float> &image, const colorscheme_params &clrpar) {
+ std::vector<uint8_t> pixels(3 * image.size());
+
+ for (size_t i = 0; i < image.size(); i++) {
+ const double iters = image[i];
+
+ std::array<uint8_t, 3> clrs = colorscheme(iters, clrpar);
+ pixels[3 * i + 0] = clrs[0];
+ pixels[3 * i + 1] = clrs[1];
+ pixels[3 * i + 2] = clrs[2];
+ }
+
+ return pixels;
+}
+
+struct job {
+ size_t from_idx, to_idx;
+};
+
+struct job_list {
+ std::vector<job> jobs;
+ std::atomic_size_t jobidx;
+};
+
+static void mandel_worker(const params &params, job_list &jl) {
+ while (true) {
+ size_t ji = jl.jobidx.fetch_add(1);
+ if (ji >= jl.jobs.size()) break;
+ mandel_job(params, jl.jobs[ji].from_idx, jl.jobs[ji].to_idx);
+ }
+}
+
+static void write_pixels_to_file(const std::string &output_fname, const std::vector<uint8_t> &pixels, size_t width, size_t height) {
+ int ret = lodepng::encode(output_fname, pixels, width, height, LCT_RGB, 8);
+ if (ret != 0) {
+ std::cerr << "Cannot write PNG: " << lodepng_error_text(ret) << std::endl;
+ exit(1);
+ }
+}
+
+int main(int argc, char **argv) {
+ size_t width = 1920, height = 1080;
+ size_t maxiter = 256;
+ window win{-2.9, 1.9, -1.35, 1.35};
+ std::string save_fname, load_fname, output_fname;
+ colorscheme_params clrpar;
+
+ for (int i = 1; i < argc; i++) {
+ const auto take_argument = [&i, argc, argv](auto &&parsefunc, auto &dest) {
+ if (i == argc - 1) {
+ std::cerr << "Flag expects argument: '" << argv[i] << "'" << std::endl;
+ exit(1);
+ }
+
+ dest = parsefunc(argv[i + 1]);
+ i++;
+ };
+
+ const auto &parse_int = [](const char *s) { return std::stoull(s); };
+ const auto &parse_flt = [](const char *s) { return std::stod(s); };
+ const auto &identity = [](const char *s) { return s; };
+
+ if (strcmp(argv[i], "-W") == 0) take_argument(parse_int, width);
+ else if (strcmp(argv[i], "-H") == 0) take_argument(parse_int, height);
+ else if (strcmp(argv[i], "-M") == 0) take_argument(parse_int, maxiter);
+ else if (strcmp(argv[i], "-o") == 0) take_argument(identity, output_fname);
+ else if (strcmp(argv[i], "-s") == 0) take_argument(identity, save_fname);
+ else if (strcmp(argv[i], "-l") == 0) take_argument(identity, load_fname);
+ else if (strcmp(argv[i], "-L") == 0) take_argument(parse_flt, win.xmin);
+ else if (strcmp(argv[i], "-R") == 0) take_argument(parse_flt, win.xmax);
+ else if (strcmp(argv[i], "-T") == 0) take_argument(parse_flt, win.ymin);
+ else if (strcmp(argv[i], "-B") == 0) take_argument(parse_flt, win.ymax);
+ else if (strcmp(argv[i], "-X") == 0) {
+ double x; take_argument(parse_flt, x);
+ win = win.recenter_x(x);
+ }
+ else if (strcmp(argv[i], "-Y") == 0) {
+ double y; take_argument(parse_flt, y);
+ win = win.recenter_y(y);
+ }
+ else if (strcmp(argv[i], "-Z") == 0) {
+ double w; take_argument(parse_flt, w);
+ win = win.set_width(w).set_height((double)height / width * w);
+ }
+ else if (strcmp(argv[i], "-cS") == 0) take_argument(parse_flt, clrpar.speed);
+ else if (strcmp(argv[i], "-cO") == 0) take_argument(parse_flt, clrpar.offset);
+ else if (strcmp(argv[i], "-h") == 0 || strcmp(argv[i], "--help") == 0) {
+ std::cout <<
+ "Usage: " << argv[0] << " [options] -o <out.png>\n"
+ " " << argv[0] << " [options] -s <data.txt>\n"
+ " " << argv[0] << " [options] -l <data.txt> -o <out.png>\n"
+ "\n"
+ " -h/--help Show this help\n"
+ "\n"
+ " -W <width> Width of output image in pixels [1920]\n"
+ " -H <height> Height of output image in pixels [1080]\n"
+ " -M <maxiter> Number of iterations after which a point is \"in set\" [256]\n"
+ " -o <out.png> Output file name\n"
+ " -s <data.txt> Iteration cache output file name (colorscheme options unused)\n"
+ " -l <data.txt> Only colorise an iteration cache file; requires -o (and\n"
+ " ignores all other flags except colorscheme flags)\n"
+ "\n"
+ "For the complex boundaries, either set them explicitly (they will potentially\n"
+ "be expanded so that the width/height ratio is satisfied):\n"
+ " -L <left> Left real bound [-2.9]\n"
+ " -R <right> Right real bound [1.9]\n"
+ " -B <bottom> Bottom imaginary bound [-1.35]\n"
+ " -T <top> Top imaginary bound [1.35]\n"
+ "or via a center point and zoom level:\n"
+ " -X <midx> Real part of midpoint of image [-0.5]\n"
+ " -Y <midy> Imaginary part of midpoint of image [0.0]\n"
+ " -Z <cplxwid> Complex width of the image (right - left) [4.8]\n"
+ "When mixing these two approaches, later arguments override earlier ones.\n"
+ "\n"
+ "Colorscheme options: (unused for -s)\n"
+ " -cS <speed> Hue spin speed [1.0]\n"
+ " -cO <offset> Hue spin offset [0.0]\n";
+ std::cout << std::flush;
+ return 0;
+ } else {
+ std::cerr << "Unrecognised option '" << argv[i] << "'" << std::endl;
+ return 1;
+ }
+ }
+
+ if ((!save_fname.empty() && (!load_fname.empty() || !output_fname.empty())) ||
+ (save_fname.empty() && output_fname.empty())) {
+ std::cerr << "Either -o, -s or -l -o is required. Use -h for help." << std::endl;
+ return 1;
+ }
+
+ if (!load_fname.empty()) {
+ std::ifstream f(load_fname, std::ios::binary);
+ if (!f) {
+ std::cerr << "Cannot open file '" << load_fname << "'" << std::endl;
+ return 1;
+ }
+ f.read((char*)&width, sizeof width);
+ f.read((char*)&height, sizeof height);
+
+ std::vector<float> image(width * height);
+ for (size_t i = 0; i < width * height; i++) {
+ f.read((char*)&image[i], sizeof(float));
+ }
+
+ std::vector<uint8_t> pixels = colorise_image(image, clrpar);
+ write_pixels_to_file(output_fname, pixels, width, height);
+ return 0;
+ }
+
+ if (width == 0 || height == 0 || maxiter == 0 ||
+ win.xmin > win.xmax || win.ymin > win.ymax) {
+ std::cerr << "Invalid parameters given" << std::endl;
+ return 1;
+ }
+
+ win = win.enlarge_to_ratio((double)width / height);
+
+ FLT xmin, xmax, ymin, ymax;
+ xmin.set_d(win.xmin);
+ xmax.set_d(win.xmax);
+ ymin.set_d(win.ymin);
+ ymax.set_d(win.ymax);
+
+ std::vector<float> image(width * height);
+
+ params params{
+ &image,
+ width, height,
+ maxiter,
+ xmin, xmax, ymin, ymax
+ };
+
+#ifdef MULTITHREAD
+ const size_t nthreads = std::thread::hardware_concurrency();
+ const size_t njobs = 16 * nthreads;
+
+ job_list jl;
+ jl.jobs.resize(njobs);
+ for (size_t i = 0; i < njobs; i++) {
+ jl.jobs[i].from_idx = i * width * height / njobs;
+ jl.jobs[i].to_idx = (i + 1) * width * height / njobs;
+ }
+ jl.jobidx.store(0);
+
+ std::vector<std::thread> threads;
+ for (size_t i = 0; i < nthreads; i++) {
+ threads.emplace_back([&params, &jl, i]() {
+ mandel_worker(params, jl);
+ });
+ }
+
+ for (std::thread &th : threads) th.join();
+#else
+ mandel_job(params, 0, width * height);
+#endif
+
+ if (!save_fname.empty()) {
+ std::ofstream f(save_fname, std::ios::binary);
+ if (!f) {
+ std::cerr << "Cannot open file '" << save_fname << "'" << std::endl;
+ return 1;
+ }
+
+ f.write((const char*)&width, sizeof width);
+ f.write((const char*)&height, sizeof height);
+ for (size_t i = 0; i < width * height; i++) {
+ f.write((const char*)&image[i], sizeof(float));
+ }
+ } else {
+ std::vector<uint8_t> pixels = colorise_image(image, clrpar);
+ write_pixels_to_file(output_fname, pixels, width, height);
+ }
+}