|
| 1 | +// The code bellow uses C++-17 features, compile it with C++-17 flags, e.g.: |
| 2 | +// clang++ -Wall -Wextra -Wshadow -Wnon-virtual-dtor -Wold-style-cast -Wcast-align -Wunused -Woverloaded-virtual -Wpedantic -Wconversion -Wsign-conversion -Wnull-dereference -Wdouble-promotion -Wformat=2 -gdwarf-3 -D_GLIBCXX_DEBUG -std=c++17 -O3 -c ./barnsley.cpp barnsley |
| 3 | + |
| 4 | +#include <array> |
| 5 | +#include <cassert> |
| 6 | +#include <fstream> |
| 7 | +#include <random> |
| 8 | + |
| 9 | +using Vec2 = std::array<double, 2>; |
| 10 | +using Vec3 = std::array<double, 3>; |
| 11 | +using Row = std::array<double, 3>; |
| 12 | +using Op = std::array<Row, 3>; |
| 13 | + |
| 14 | +constexpr auto OpN = 4U; |
| 15 | + |
| 16 | +template <size_t N> |
| 17 | +auto operator+(std::array<double, N> x, std::array<double, N> y) { |
| 18 | + for (auto i = 0U; i < N; ++i) |
| 19 | + x[i] += y[i]; |
| 20 | + return x; |
| 21 | +} |
| 22 | + |
| 23 | +template <size_t N> |
| 24 | +auto operator*(double k, std::array<double, N> v) { |
| 25 | + for (auto i = 0U; i < N; ++i) |
| 26 | + v[i] *= k; |
| 27 | + return v; |
| 28 | +} |
| 29 | + |
| 30 | +template <size_t N> |
| 31 | +auto operator*(std::array<double, N> v, double k) { |
| 32 | + return k * v; |
| 33 | +} |
| 34 | + |
| 35 | +auto operator*(const Op& x, const Vec3& y) { |
| 36 | + auto ret = Vec3{}; |
| 37 | + for (auto i = 0U; i < 3U; ++i) { |
| 38 | + ret[i] = 0; |
| 39 | + for (auto j = 0U; j < 3U; ++j) |
| 40 | + ret[i] += y[j] * x[i][j]; |
| 41 | + } |
| 42 | + return ret; |
| 43 | +} |
| 44 | + |
| 45 | +// Returns a pseudo-random number generator |
| 46 | +std::default_random_engine& rng() { |
| 47 | + // Initialize static pseudo-random engine with non-deterministic random seed |
| 48 | + static std::default_random_engine randEngine(std::random_device{}()); |
| 49 | + return randEngine; |
| 50 | +} |
| 51 | + |
| 52 | +// Returns a random double in [0, 1) |
| 53 | +double drand() { |
| 54 | + return std::uniform_real_distribution<double>(0.0, 1.0)(rng()); |
| 55 | +} |
| 56 | + |
| 57 | +// This is a function that reads in the Hutchinson operator and |
| 58 | +// corresponding |
| 59 | +// probabilities and outputs a randomly selected transform |
| 60 | +// This works by choosing a random number and then iterating through all |
| 61 | +// probabilities until it finds an appropriate bin |
| 62 | +auto select_array( |
| 63 | + const std::array<Op, OpN>& hutchinson_op, |
| 64 | + const std::array<double, OpN>& probabilities) { |
| 65 | + |
| 66 | + // random number to be binned |
| 67 | + auto rnd = drand(); |
| 68 | + |
| 69 | + // This checks to see if a random number is in a bin, if not, that |
| 70 | + // probability is subtracted from the random number and we check the |
| 71 | + // next bin in the list |
| 72 | + for (auto i = 0U; i < probabilities.size(); ++i) { |
| 73 | + if (rnd < probabilities[i]) |
| 74 | + return hutchinson_op[i]; |
| 75 | + rnd -= probabilities[i]; |
| 76 | + } |
| 77 | + assert(!static_cast<bool>("check if probabilities adding up to 1")); |
| 78 | +} |
| 79 | + |
| 80 | +// This is a general function to simulate a chaos game |
| 81 | +// n is the number of iterations |
| 82 | +// initial_location is the the starting point of the chaos game |
| 83 | +// hutchinson_op is the set of functions to iterate through |
| 84 | +// probabilities is the set of probabilities corresponding to the likelihood |
| 85 | +// of choosing their corresponding function in hutchinson_op |
| 86 | +auto chaos_game( |
| 87 | + size_t n, |
| 88 | + Vec2 initial_location, |
| 89 | + const std::array<Op, OpN>& hutchinson_op, |
| 90 | + const std::array<double, OpN>& probabilities) { |
| 91 | + |
| 92 | + // Initializing the output array and the initial point |
| 93 | + auto output_points = std::vector<Vec2>{}; |
| 94 | + |
| 95 | + // extending point to 3D for affine transform |
| 96 | + auto point = Vec3{initial_location[0], initial_location[1], 1}; |
| 97 | + |
| 98 | + for (auto i = 0U; i < n; ++i) { |
| 99 | + output_points.push_back(Vec2{point[0], point[1]}); |
| 100 | + point = select_array(hutchinson_op, probabilities) * point; |
| 101 | + } |
| 102 | + |
| 103 | + return output_points; |
| 104 | +} |
| 105 | + |
| 106 | +int main() { |
| 107 | + |
| 108 | + const std::array barnsley_hutchinson = { |
| 109 | + Op{Row{0.0, 0.0, 0.0}, Row{0.0, 0.16, 0.0}, Row{0.0, 0.0, 1.0}}, |
| 110 | + Op{Row{0.85, 0.04, 0.0}, Row{-0.04, 0.85, 1.60}, Row{0.0, 0.0, 1.0}}, |
| 111 | + Op{Row{0.20, -0.26, 0.0}, Row{0.23, 0.22, 1.60}, Row{0.0, 0.0, 1.0}}, |
| 112 | + Op{Row{-0.15, 0.28, 0.0}, Row{0.26, 0.24, 0.44}, Row{0.0, 0.0, 1.0}}}; |
| 113 | + |
| 114 | + const std::array barnsley_probabilities = {0.01, 0.85, 0.07, 0.07}; |
| 115 | + auto output_points = chaos_game( |
| 116 | + 10'000, Vec2{0, 0}, barnsley_hutchinson, barnsley_probabilities); |
| 117 | + |
| 118 | + std::ofstream ofs("out.dat"); |
| 119 | + for (auto pt : output_points) |
| 120 | + ofs << pt[0] << '\t' << pt[1] << '\n'; |
| 121 | +} |
0 commit comments